In [2]:

    

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


    

In [3]:

    

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 [4]:

    

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


    

In [6]:

    

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


    

    
Out[6]:





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

In [7]:

    

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


    

    
Out[7]:





$$- c^{2} p - c \gamma + p^{5} v_{p} + p^{3} + p$$

In [8]:

    

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


    

    
Out[8]:





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

In [9]:

    

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


    

    
Out[9]:





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

In [10]:

    

series(E_c,c,n=7)


    

    
Out[10]:





$$c \left(- \gamma + \frac{r_{c}}{\gamma}\right) + c^{3} \left(- \frac{r_{c}}{\gamma} + \frac{u_{c}}{\gamma} + \frac{r_{c}^{3}}{\gamma^{3}} - \frac{r_{c}^{2}}{\gamma^{3}}\right) + c^{5} \left(- \frac{u_{c}}{\gamma} + \frac{v_{c}}{\gamma} + \frac{3 u_{c}}{\gamma^{3}} r_{c}^{2} + \frac{r_{c}^{2}}{\gamma^{3}} - \frac{2 r_{c}}{\gamma^{3}} u_{c} + \frac{r_{c}^{5} v_{p}}{\gamma^{5}} - \frac{3 r_{c}^{4}}{\gamma^{5}} + \frac{2 r_{c}^{3}}{\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 [13]:

    

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

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

Q = (-u_c/gamma+v_c/gamma+3*u_c*r_c**2/gamma**3+r_c**2/gamma**3
     -2*r_c*u_c/gamma**3+v_p*(r_c/gamma)**5-3*r_c**4/gamma**5+2*r_c**3/gamma**5)


    

In [19]:

    

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


    

    
Out[19]:





$$a^{3} + a^{2} \left(3 \gamma - \frac{1}{\gamma}\right) + a \left(3 \gamma^{2} - 3\right) + \gamma^{3} - 2 \gamma + \frac{u_{c}}{\gamma}$$

In [20]:

    

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


    

    
Out[20]:





$$a^{5} v_{p} + a^{4} \left(5 \gamma v_{p} - \frac{3}{\gamma}\right) + a^{3} \left(10 \gamma^{2} v_{p} - 12 + \frac{2}{\gamma^{2}}\right) + a^{2} \left(10 \gamma^{3} v_{p} - 18 \gamma + \frac{3 u_{c}}{\gamma} + \frac{7}{\gamma}\right) + a \left(5 \gamma^{4} v_{p} - 12 \gamma^{2} + 6 u_{c} + 8 - \frac{2 u_{c}}{\gamma^{2}}\right) + \gamma^{5} v_{p} - 3 \gamma^{3} + 3 \gamma u_{c} + 3 \gamma - \frac{3 u_{c}}{\gamma} + \frac{v_{c}}{\gamma}$$

In [18]:

    

b0 = gamma**3-2*gamma+u_c/gamma

b1 = 3*gamma**2-3

b2 = 3*gamma-1/gamma

b3 = 1

q0 = gamma**5*v_p-3*gamma**3+3*gamma*u_c+3*gamma-3*u_c/gamma+v_c/gamma

q1 = 5*gamma**4*v_p-12*gamma**2+6*u_c+8-2*u_c/gamma**2

q2 = 10*v_p*gamma**3-18*gamma+(7+3*u_c)/gamma

q3 = 10*v_p*gamma**2-12+2/gamma**2

q4 = 5*v_p*gamma-3/gamma

q5 = v_p


    

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 [ ]: