In [1]:

    

%matplotlib inline
from sympy import *
init_printing(use_unicode=True)


    

In [2]:

    

r, u, v, c, r_c, u_c, v_c, E, p, r_p, u_p, v_p, e, a, b, q, b_0, b_1, b_2, b_3, q_0, q_1, q_2, q_3, q_4, q_5  = symbols('r u v c r_c u_c v_c E p r_p u_p v_p e a b q b_0 b_1 b_2 b_3 q_0 q_1 q_2 q_3 q_4 q_5')


    

In [3]:

    

gamma = symbols('gamma',positive=True)


    

In [4]:

    

f2 = ((1/2)*r_c*c**2+(1/4)*u_c*c**4+(1/6)*v_c*c**6+(1/2)*r_p*p**2
      +(1/4)*u_p*p**4+(1/6)*v_p*p**6-E*p-gamma*c*p-e*c**2*p**2)
nsimplify(f2)


    

    
Out[4]:





$$- E p + \frac{c^{6} v_{c}}{6} + \frac{c^{4} u_{c}}{4} - c^{2} e p^{2} + \frac{c^{2} r_{c}}{2} - c \gamma p + \frac{p^{6} v_{p}}{6} + \frac{p^{4} u_{p}}{4} + \frac{p^{2} r_{p}}{2}$$

In [5]:

    

E_c = solve(f2.diff(p),E)[0]
E_c


    

    
Out[5]:





$$- 2.0 c^{2} e p - c \gamma + p^{5} v_{p} + p^{3} u_{p} + p r_{p}$$

In [6]:

    

p_min = nsimplify(solve(f2.diff(c),p)[0])
p_min


    

    
Out[6]:





$$\frac{1}{4 c e} \left(- \gamma + \sqrt{8 c^{6} e v_{c} + 8 c^{4} e u_{c} + 8 c^{2} e r_{c} + \gamma^{2}}\right)$$

In [7]:

    

E_c = nsimplify(E_c.subs(p,p_min))
E_c


    

    
Out[7]:





$$- c \gamma - \frac{c}{2} \left(- \gamma + \sqrt{8 c^{6} e v_{c} + 8 c^{4} e u_{c} + 8 c^{2} e r_{c} + \gamma^{2}}\right) + \frac{r_{p}}{4 c e} \left(- \gamma + \sqrt{8 c^{6} e v_{c} + 8 c^{4} e u_{c} + 8 c^{2} e r_{c} + \gamma^{2}}\right) + \frac{u_{p}}{64 c^{3} e^{3}} \left(- \gamma + \sqrt{8 c^{6} e v_{c} + 8 c^{4} e u_{c} + 8 c^{2} e r_{c} + \gamma^{2}}\right)^{3} + \frac{v_{p}}{1024 c^{5} e^{5}} \left(- \gamma + \sqrt{8 c^{6} e v_{c} + 8 c^{4} e u_{c} + 8 c^{2} e r_{c} + \gamma^{2}}\right)^{5}$$

In [8]:

    

series(E_c,c,n=7)


    

    
Out[8]:





$$c \left(- \gamma + \frac{r_{c} r_{p}}{\gamma}\right) + c^{3} \left(- \frac{2 e}{\gamma} r_{c} - \frac{2 e}{\gamma^{3}} r_{c}^{2} r_{p} + \frac{r_{p} u_{c}}{\gamma} + \frac{r_{c}^{3} u_{p}}{\gamma^{3}}\right) + c^{5} \left(\frac{4 e^{2}}{\gamma^{3}} r_{c}^{2} + \frac{8 r_{p}}{\gamma^{5}} e^{2} r_{c}^{3} - \frac{2 e}{\gamma} u_{c} - \frac{4 e}{\gamma^{3}} r_{c} r_{p} u_{c} - \frac{6 e}{\gamma^{5}} r_{c}^{4} u_{p} + \frac{r_{p} v_{c}}{\gamma} + \frac{3 u_{c}}{\gamma^{3}} r_{c}^{2} u_{p} + \frac{r_{c}^{5} v_{p}}{\gamma^{5}}\right) + \mathcal{O}\left(c^{7}\right)$$

In [9]:

    

Etrun = a*c+b*c**3+q*c**5


    

In [10]:

    

solve(Etrun.diff(c),c)


    

    
Out[10]:





$$\left [ - \frac{\sqrt{10}}{10} \sqrt{- \frac{3 b}{q} - \frac{1}{q} \sqrt{- 20 a q + 9 b^{2}}}, \quad \frac{\sqrt{10}}{10} \sqrt{- \frac{3 b}{q} - \frac{1}{q} \sqrt{- 20 a q + 9 b^{2}}}, \quad - \frac{\sqrt{10}}{10} \sqrt{- \frac{3 b}{q} + \frac{1}{q} \sqrt{- 20 a q + 9 b^{2}}}, \quad \frac{\sqrt{10}}{10} \sqrt{- \frac{3 b}{q} + \frac{1}{q} \sqrt{- 20 a q + 9 b^{2}}}\right ]$$

In [11]:

    

c_L = solve(Etrun.diff(c),c)[1]
c_U = solve(Etrun.diff(c),c)[3]
c_L,c_U


    

    
Out[11]:





$$\left ( \frac{\sqrt{10}}{10} \sqrt{- \frac{3 b}{q} - \frac{1}{q} \sqrt{- 20 a q + 9 b^{2}}}, \quad \frac{\sqrt{10}}{10} \sqrt{- \frac{3 b}{q} + \frac{1}{q} \sqrt{- 20 a q + 9 b^{2}}}\right )$$

In [12]:

    

E_L = simplify(Etrun.subs(c,c_U))
E_U = simplify(Etrun.subs(c,c_L))


    

In [13]:

    

E_L,E_U


    

    
Out[13]:





$$\left ( \frac{\sqrt{10}}{250 q} \sqrt{\frac{1}{q} \left(- 3 b + \sqrt{- 20 a q + 9 b^{2}}\right)} \left(20 a q - 3 b^{2} + b \sqrt{- 20 a q + 9 b^{2}}\right), \quad \frac{\sqrt{10}}{250 q} \sqrt{\frac{1}{q} \left(- 3 b - \sqrt{- 20 a q + 9 b^{2}}\right)} \left(20 a q - 3 b^{2} - b \sqrt{- 20 a q + 9 b^{2}}\right)\right )$$

In [8]:

    

rc = (gamma**2+a*gamma)/r_p

B = (r_p*u_c/gamma+u_p*(r_c/gamma)**3-2*e*r_c/gamma-2*e*r_p*r_c**2/gamma**3)

Q = ((2*e*r_c)**2/gamma**3+8*r_p*e**2*r_c**3/gamma**5-2*e*u_c/gamma
     -4*e*r_c*r_p*u_c/gamma**3-6*e*u_p*r_c**4/gamma**5+r_p*v_c/gamma
     +3*u_c*u_p*r_c**2/gamma**3+v_p*(r_c/gamma)**5)


    

In [9]:

    

collect(expand(B.subs(r_c,rc)),a)


    

    
Out[9]:





$$\frac{a^{3} u_{p}}{r_{p}^{3}} + a^{2} \left(- \frac{2 e}{\gamma r_{p}} + \frac{3 u_{p}}{r_{p}^{3}} \gamma\right) + a \left(- \frac{6 e}{r_{p}} + \frac{3 u_{p}}{r_{p}^{3}} \gamma^{2}\right) - \frac{4 e}{r_{p}} \gamma + \frac{\gamma^{3} u_{p}}{r_{p}^{3}} + \frac{r_{p} u_{c}}{\gamma}$$

In [10]:

    

collect(expand(Q.subs(r_c,rc)),a)


    

    
Out[10]:





$$\frac{a^{5} v_{p}}{r_{p}^{5}} + a^{4} \left(- \frac{6 e u_{p}}{\gamma r_{p}^{4}} + \frac{5 v_{p}}{r_{p}^{5}} \gamma\right) + a^{3} \left(\frac{8 e^{2}}{\gamma^{2} r_{p}^{2}} - \frac{24 e}{r_{p}^{4}} u_{p} + \frac{10 v_{p}}{r_{p}^{5}} \gamma^{2}\right) + a^{2} \left(\frac{28 e^{2}}{\gamma r_{p}^{2}} - \frac{36 e}{r_{p}^{4}} \gamma u_{p} + \frac{10 v_{p}}{r_{p}^{5}} \gamma^{3} + \frac{3 u_{c} u_{p}}{\gamma r_{p}^{2}}\right) + a \left(\frac{32 e^{2}}{r_{p}^{2}} - \frac{24 e}{r_{p}^{4}} \gamma^{2} u_{p} - \frac{4 e}{\gamma^{2}} u_{c} + \frac{5 v_{p}}{r_{p}^{5}} \gamma^{4} + \frac{6 u_{c}}{r_{p}^{2}} u_{p}\right) + \frac{12 \gamma}{r_{p}^{2}} e^{2} - \frac{6 e}{r_{p}^{4}} \gamma^{3} u_{p} - \frac{6 e}{\gamma} u_{c} + \frac{\gamma^{5} v_{p}}{r_{p}^{5}} + \frac{3 u_{c}}{r_{p}^{2}} \gamma u_{p} + \frac{r_{p} v_{c}}{\gamma}$$

In [11]:

    

b0 = (u_p*(gamma/r_p)**3-4*e*gamma/r_p+u_c*r_p/gamma)

b1 = (3*u_p*gamma**2/r_p**3-6*e/r_p)

b2 = (3*u_p*gamma/r_p**3-2*e/(gamma*r_p))

b3 = u_p/r_p**3

q0 = (v_p*(gamma/r_p)**5-6*e*u_p*gamma**3/r_p**4+3*(gamma/r_p**2)*(4*e**2+u_c*u_p)-6*e*u_c/gamma+v_c*r_p/gamma)

q1 = (5*v_p*gamma**4/r_p**5-24*e*u_p*gamma**2/r_p**4+2*(16*e**2+3*u_c*u_p)/r_p**2-4*e*u_c/gamma**2)

q2 = (10*v_p*gamma**3/r_p**5-36*e*u_p*gamma/r_p**4+(28*e**2+3*u_c*u_p)/(gamma*r_p**2))

q3 = (10*v_p*gamma**2/r_p**5-24*e*u_p/r_p**4+8*(e/(gamma*r_p))**2)

q4 = (5*v_p*gamma/r_p**5-6*e*u_p/(r_p**4*gamma))

q5 = v_p/r_p**5


    

In [12]:

    

uc = solve(b0-b_0,u_c)[0]
up = solve(b3-b_3,u_p)[0]
vc = solve(q0-q_0,v_c)[0]
vp = solve(q5-q_5,v_p)[0]


    

In [13]:

    

replacements = [(v_p,vp),(u_p,up),(u_c,uc)]


    

In [16]:

    

vc = simplify(vc.subs([i for i in replacements]))


    

In [17]:

    

expand(vc)


    

    
Out[17]:





$$- \frac{3 b_{0}}{r_{p}} b_{3} \gamma^{3} + \frac{6 b_{0}}{r_{p}^{2}} e \gamma + \frac{3 b_{3}^{2}}{r_{p}} \gamma^{6} - \frac{12 b_{3}}{r_{p}^{2}} e \gamma^{4} + \frac{12 e^{2}}{r_{p}^{3}} \gamma^{2} - \frac{\gamma^{6} q_{5}}{r_{p}} + \frac{\gamma q_{0}}{r_{p}}$$

In [18]:

    

uc = simplify(uc.subs([i for i in replacements]))


    

In [19]:

    

expand(uc)


    

    
Out[19]:





$$\frac{b_{0} \gamma}{r_{p}} - \frac{b_{3} \gamma^{4}}{r_{p}} + \frac{4 e}{r_{p}^{2}} \gamma^{2}$$

In [20]:

    

up


    

    
Out[20]:





$$b_{3} r_{p}^{3}$$

In [21]:

    

vp


    

    
Out[21]:





$$q_{5} r_{p}^{5}$$

In [22]:

    

b0


    

    
Out[22]:





$$- \frac{4 e}{r_{p}} \gamma + \frac{\gamma^{3} u_{p}}{r_{p}^{3}} + \frac{r_{p} u_{c}}{\gamma}$$

In [23]:

    

b1 = simplify(b1.subs(u_p,up))


    

In [24]:

    

b1


    

    
Out[24]:





$$3 b_{3} \gamma^{2} - \frac{6 e}{r_{p}}$$

In [25]:

    

b2 = simplify(b2.subs(u_p,up))


    

In [26]:

    

b2


    

    
Out[26]:





$$3 b_{3} \gamma - \frac{2 e}{\gamma r_{p}}$$

In [27]:

    

b3


    

    
Out[27]:





$$\frac{u_{p}}{r_{p}^{3}}$$

In [28]:

    

B_a = b_0+b1*a+b2*a**2+b_3*a**3


    

In [29]:

    

B_a


    

    
Out[29]:





$$a^{3} b_{3} + a^{2} \left(3 b_{3} \gamma - \frac{2 e}{\gamma r_{p}}\right) + a \left(3 b_{3} \gamma^{2} - \frac{6 e}{r_{p}}\right) + b_{0}$$

In [30]:

    

q0


    

    
Out[30]:





$$- \frac{6 e}{r_{p}^{4}} \gamma^{3} u_{p} - \frac{6 e}{\gamma} u_{c} + \frac{\gamma^{5} v_{p}}{r_{p}^{5}} + \frac{3 \gamma}{r_{p}^{2}} \left(4 e^{2} + u_{c} u_{p}\right) + \frac{r_{p} v_{c}}{\gamma}$$

In [33]:

    

q1 = simplify(q1.subs([i for i in replacements]))


    

In [34]:

    

q1


    

    
Out[34]:





$$6 b_{0} b_{3} \gamma - \frac{4 b_{0} e}{\gamma r_{p}} - 6 b_{3}^{2} \gamma^{4} + \frac{4 b_{3}}{r_{p}} e \gamma^{2} + \frac{16 e^{2}}{r_{p}^{2}} + 5 \gamma^{4} q_{5}$$

In [37]:

    

q2 = q2.subs([i for i in replacements])


    

In [38]:

    

expand(q2)


    

    
Out[38]:





$$3 b_{0} b_{3} - 3 b_{3}^{2} \gamma^{3} - \frac{24 b_{3}}{r_{p}} e \gamma + \frac{28 e^{2}}{\gamma r_{p}^{2}} + 10 \gamma^{3} q_{5}$$

In [39]:

    

q3 = simplify(q3.subs([i for i in replacements]))


    

In [40]:

    

q3


    

    
Out[40]:





$$- \frac{24 b_{3}}{r_{p}} e + \frac{8 e^{2}}{\gamma^{2} r_{p}^{2}} + 10 \gamma^{2} q_{5}$$

In [41]:

    

q4 = simplify(q4.subs([i for i in replacements]))


    

In [42]:

    

q4


    

    
Out[42]:





$$- \frac{6 b_{3} e}{\gamma r_{p}} + 5 \gamma q_{5}$$

In [43]:

    

Q_a = q_0+q1*a+q2*a**2+q3*a**3+q4*a**4+q_5*a**5


    

In [44]:

    

collect(expand(Q_a),a)


    

    
Out[44]:





$$a^{5} q_{5} + a^{4} \left(- \frac{6 b_{3} e}{\gamma r_{p}} + 5 \gamma q_{5}\right) + a^{3} \left(- \frac{24 b_{3}}{r_{p}} e + \frac{8 e^{2}}{\gamma^{2} r_{p}^{2}} + 10 \gamma^{2} q_{5}\right) + a^{2} \left(3 b_{0} b_{3} - 3 b_{3}^{2} \gamma^{3} - \frac{24 b_{3}}{r_{p}} e \gamma + \frac{28 e^{2}}{\gamma r_{p}^{2}} + 10 \gamma^{3} q_{5}\right) + a \left(6 b_{0} b_{3} \gamma - \frac{4 b_{0} e}{\gamma r_{p}} - 6 b_{3}^{2} \gamma^{4} + \frac{4 b_{3}}{r_{p}} e \gamma^{2} + \frac{16 e^{2}}{r_{p}^{2}} + 5 \gamma^{4} q_{5}\right) + q_{0}$$

In [45]:

    

series(B_a**2-(20*a*Q_a/9),a,n=7)


    

    
Out[45]:





$$a^{6} \left(b_{3}^{2} - \frac{20 q_{5}}{9}\right) + a^{5} \left(6 b_{3}^{2} \gamma + \frac{28 b_{3} e}{3 \gamma r_{p}} - \frac{100 q_{5}}{9} \gamma\right) + a^{4} \left(15 b_{3}^{2} \gamma^{2} + \frac{88 b_{3} e}{3 r_{p}} - \frac{124 e^{2}}{9 \gamma^{2} r_{p}^{2}} - \frac{200 q_{5}}{9} \gamma^{2}\right) + a^{3} \left(- \frac{14 b_{0}}{3} b_{3} + \frac{74 b_{3}^{2}}{3} \gamma^{3} + \frac{16 b_{3} e \gamma}{3 r_{p}} - \frac{344 e^{2}}{9 \gamma r_{p}^{2}} - \frac{200 q_{5}}{9} \gamma^{3}\right) + a^{2} \left(- \frac{22 b_{0}}{3} b_{3} \gamma + \frac{44 b_{0} e}{9 \gamma r_{p}} + \frac{67 b_{3}^{2}}{3} \gamma^{4} - \frac{404 b_{3} e \gamma^{2}}{9 r_{p}} + \frac{4 e^{2}}{9 r_{p}^{2}} - \frac{100 q_{5}}{9} \gamma^{4}\right) + a \left(6 b_{0} b_{3} \gamma^{2} - \frac{12 b_{0}}{r_{p}} e - \frac{20 q_{0}}{9}\right) + b_{0}^{2}$$

In [ ]: