In [2]:

    

from sympy import *
init_session()


    

    




IPython console for SymPy 0.7.4.1-git (Python 3.3.3-64-bit) (ground types: python)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)

Documentation can be found at http://www.sympy.org

In [3]:

    

alpha, beta, sigma, epsilon, deltaA, deltaM, deltaE = symbols('alpha beta sigma epsilon delta_A delta_M delta_E')
k1, k2 = symbols('k1 k2')
a, R, deltaP = symbols('a R delta_P')
B = 1


    

In [4]:

    

A, M, E = Function("A")(t), Function("M")(t), Function("E")(t)


    

In [5]:

    

f1 = lambda p: p**3/(k1 + p**3)
f2 = lambda p: a*k2**2/(k2**2 + p**2)
p = (R*B*E)/deltaP


    

In [6]:

    

dA = (sigma + alpha*M)*f1(p) - (beta + deltaA)*A - epsilon*A**2
dM = 8*beta*f2(p)*A - f1(p)*alpha*M - deltaM*M 
dE = 60*beta*(1 - f2(p))*A - deltaE*E


    

In [7]:

    

dA


    

    
Out[7]:





$$\frac{R^{3} \left(\alpha M{\left (t \right )} + \sigma\right) E^{3}{\left (t \right )}}{\delta_{P}^{3} \left(\frac{R^{3}}{\delta_{P}^{3}} E^{3}{\left (t \right )} + k_{1}\right)} - \epsilon A^{2}{\left (t \right )} - \left(\beta + \delta_{A}\right) A{\left (t \right )}$$

In [8]:

    

dM


    

    
Out[8]:





$$- \frac{R^{3} \alpha E^{3}{\left (t \right )} M{\left (t \right )}}{\delta_{P}^{3} \left(\frac{R^{3}}{\delta_{P}^{3}} E^{3}{\left (t \right )} + k_{1}\right)} + \frac{8 a \beta k_{2}^{2} A{\left (t \right )}}{\frac{R^{2}}{\delta_{P}^{2}} E^{2}{\left (t \right )} + k_{2}^{2}} - \delta_{M} M{\left (t \right )}$$

In [9]:

    

dE


    

    
Out[9]:





$$60 \beta \left(- \frac{a k_{2}^{2}}{\frac{R^{2}}{\delta_{P}^{2}} E^{2}{\left (t \right )} + k_{2}^{2}} + 1\right) A{\left (t \right )} - \delta_{E} E{\left (t \right )}$$

In [19]:

    

A, M, E = symbols('A M E')
k1 = 2
k2 = 1
m = 3
a = 0.7
sigma = 0.02
alpha = 20
epsilon = 1
deltaM = 0.01
deltaP = 1
deltaE = 0.3
R = 50
f1 = lambda p: p**3/(k1 + p**3)
f2 = lambda p: a*k2**2/(k2**2 + p**2)
p = (R*B*E)/deltaP
lm = Matrix([(sigma + alpha*M)*f1(p) - (1)*A - epsilon*A**2, f2(p)*A - f1(p)*alpha*M - deltaM*M , 0.1*(1 - f2(p))*A - deltaE*E])
J = lm.jacobian(Matrix([A, M, E]))


    

In [25]:

    

[N(temp) for temp in J.subs({ A: 0.07016951, M: 0.007120757, E: 0.02117123 }).eigenvals().keys()]


    

    
Out[25]:





$$\begin{bmatrix}-0.625038327883738 + 1.0 \cdot 10^{-22} i, & -0.275557641844199 - 1.0 \cdot 10^{-22} i, & -7.87987121624758 - 7.0 \cdot 10^{-24} i\end{bmatrix}$$

In [45]:

    

charpoly = lm.berkowitz_charpoly().as_expr()
charpoly


    

    
Out[45]:





$$\lambda^{3} + \frac{\lambda^{2}}{62500 E^{3}{\left (t \right )} + 1} \left(\left(A{\left (t \right )} + 1.31\right) \left(62500 E^{3}{\left (t \right )} + 1\right) + 1250000 E^{3}{\left (t \right )}\right) + \frac{1.0 \lambda}{312500000.0 E^{5}{\left (t \right )} + 125000.0 E^{3}{\left (t \right )} + 5000.0 E^{2}{\left (t \right )} + 2.0} \left(6346875000.0 A{\left (t \right )} E^{5}{\left (t \right )} + 2538750.0 A{\left (t \right )} E^{3}{\left (t \right )} + 1550.0 A{\left (t \right )} E^{2}{\left (t \right )} + 0.62 A{\left (t \right )} + 8222812500.0 E^{5}{\left (t \right )} - 31250000.0 E^{4}{\left (t \right )} + 1539125.0 E^{3}{\left (t \right )} - 2185.0 E^{2}{\left (t \right )} + 0.626\right) + \frac{1}{\left(2500 E^{2}{\left (t \right )} + 1\right) \left(62500 E^{3}{\left (t \right )} + 1\right)^{2}} \left(\frac{1}{2} \left(125000 \left(A{\left (t \right )} + 1\right) \left(250.0 E^{2}{\left (t \right )} + 0.03\right) - 525000.0 E{\left (t \right )}\right) \left(62500 E^{3}{\left (t \right )} + 1\right) E^{2}{\left (t \right )} - 62500 \left(\left(A{\left (t \right )} + 1.01\right) \left(62500 E^{3}{\left (t \right )} + 1\right) + 1250000 E^{3}{\left (t \right )}\right) \left(250.0 E^{2}{\left (t \right )} + 0.03\right) E^{2}{\left (t \right )} + 0.15 \left(A{\left (t \right )} + 1\right) \left(2500 E^{2}{\left (t \right )} + 1\right) \left(62500 E^{3}{\left (t \right )} + 1\right) \left(2501250.0 E^{3}{\left (t \right )} + 0.02\right)\right)$$

In [47]:

    

solve(charpoly.subs({A: 0.06167574, M: 0.00599124, E: 0.01866815}))


    

    
Out[47]:





$$\begin{bmatrix}-6.22412163098138 - 5.0 \cdot 10^{-23} i, & -1.42654569846416 + 2.0 \cdot 10^{-22} i, & 0.497518892083428 - 1.0 \cdot 10^{-22} i\end{bmatrix}$$

In [ ]:

    




    

In [66]:

    

ep = Symbol('Epsilon')
var('a:10 m:10 e:10')
Abar = a0(t) + ep*a1(t) + ep**2*a2(t)
Mbar = m0(t) + ep*m1(t) + ep**2*m2(t)
Ebar = a0(t) + ep*e1(t) + ep**2*e2(t)


    

In [67]:

    

dAbar, dMbar, dEbar = [f.subs({A: Abar, M: Mbar, E: Ebar}) for f in (dA, dM, dE)]


    

In [68]:

    

# Simplify, get rid of denominator since equations equal 0
dAbar = fraction((dAbar - Abar.diff(t)).expand().combsimp())[0].collect(ep)
dMbar = fraction((dMbar - Mbar.diff(t)).expand().combsimp())[0].collect(ep)
dEbar = fraction((dMbar - Ebar.diff(t)).expand().combsimp())[0].collect(ep)


    

In [72]:

    

dAbar.coeff(ep, 0).simplify()


    

    
Out[72]:





$$R^{3} \alpha \operatorname{a0}^{3}{\left (t \right )} \operatorname{m0}{\left (t \right )} - R^{3} \beta \operatorname{a0}^{4}{\left (t \right )} - R^{3} \delta_{A} \operatorname{a0}^{4}{\left (t \right )} - R^{3} \epsilon \operatorname{a0}^{5}{\left (t \right )} + R^{3} \sigma \operatorname{a0}^{3}{\left (t \right )} - R^{3} \operatorname{a0}^{3}{\left (t \right )} \frac{d}{d t} \operatorname{a0}{\left (t \right )} - \beta \delta_{P}^{3} k_{1} \operatorname{a0}{\left (t \right )} - \delta_{A} \delta_{P}^{3} k_{1} \operatorname{a0}{\left (t \right )} - \delta_{P}^{3} \epsilon k_{1} \operatorname{a0}^{2}{\left (t \right )} - \delta_{P}^{3} k_{1} \frac{d}{d t} \operatorname{a0}{\left (t \right )}$$

In [73]:

    

dAbar.coeff(ep, 1).simplify()


    

    
Out[73]:





$$R^{3} \alpha \operatorname{a0}^{3}{\left (t \right )} \operatorname{m1}{\left (t \right )} + 3 R^{3} \alpha \operatorname{a0}^{2}{\left (t \right )} \operatorname{e1}{\left (t \right )} \operatorname{m0}{\left (t \right )} - R^{3} \beta \operatorname{a0}^{3}{\left (t \right )} \operatorname{a1}{\left (t \right )} - 3 R^{3} \beta \operatorname{a0}^{3}{\left (t \right )} \operatorname{e1}{\left (t \right )} - R^{3} \delta_{A} \operatorname{a0}^{3}{\left (t \right )} \operatorname{a1}{\left (t \right )} - 3 R^{3} \delta_{A} \operatorname{a0}^{3}{\left (t \right )} \operatorname{e1}{\left (t \right )} - 2 R^{3} \epsilon \operatorname{a0}^{4}{\left (t \right )} \operatorname{a1}{\left (t \right )} - 3 R^{3} \epsilon \operatorname{a0}^{4}{\left (t \right )} \operatorname{e1}{\left (t \right )} + 3 R^{3} \sigma \operatorname{a0}^{2}{\left (t \right )} \operatorname{e1}{\left (t \right )} - R^{3} \operatorname{a0}^{3}{\left (t \right )} \frac{d}{d t} \operatorname{a1}{\left (t \right )} - 3 R^{3} \operatorname{a0}^{2}{\left (t \right )} \operatorname{e1}{\left (t \right )} \frac{d}{d t} \operatorname{a0}{\left (t \right )} - \beta \delta_{P}^{3} k_{1} \operatorname{a1}{\left (t \right )} - \delta_{A} \delta_{P}^{3} k_{1} \operatorname{a1}{\left (t \right )} - 2 \delta_{P}^{3} \epsilon k_{1} \operatorname{a0}{\left (t \right )} \operatorname{a1}{\left (t \right )} - \delta_{P}^{3} k_{1} \frac{d}{d t} \operatorname{a1}{\left (t \right )}$$

In [74]:

    

dAbar.coeff(ep, 2).simplify()


    

    
Out[74]:





$$R^{3} \alpha \operatorname{a0}^{3}{\left (t \right )} \operatorname{m2}{\left (t \right )} + 3 R^{3} \alpha \operatorname{a0}^{2}{\left (t \right )} \operatorname{e1}{\left (t \right )} \operatorname{m1}{\left (t \right )} + 3 R^{3} \alpha \operatorname{a0}^{2}{\left (t \right )} \operatorname{e2}{\left (t \right )} \operatorname{m0}{\left (t \right )} + 3 R^{3} \alpha \operatorname{a0}{\left (t \right )} \operatorname{e1}^{2}{\left (t \right )} \operatorname{m0}{\left (t \right )} - R^{3} \beta \operatorname{a0}^{3}{\left (t \right )} \operatorname{a2}{\left (t \right )} - 3 R^{3} \beta \operatorname{a0}^{3}{\left (t \right )} \operatorname{e2}{\left (t \right )} - 3 R^{3} \beta \operatorname{a0}^{2}{\left (t \right )} \operatorname{a1}{\left (t \right )} \operatorname{e1}{\left (t \right )} - 3 R^{3} \beta \operatorname{a0}^{2}{\left (t \right )} \operatorname{e1}^{2}{\left (t \right )} - R^{3} \delta_{A} \operatorname{a0}^{3}{\left (t \right )} \operatorname{a2}{\left (t \right )} - 3 R^{3} \delta_{A} \operatorname{a0}^{3}{\left (t \right )} \operatorname{e2}{\left (t \right )} - 3 R^{3} \delta_{A} \operatorname{a0}^{2}{\left (t \right )} \operatorname{a1}{\left (t \right )} \operatorname{e1}{\left (t \right )} - 3 R^{3} \delta_{A} \operatorname{a0}^{2}{\left (t \right )} \operatorname{e1}^{2}{\left (t \right )} - 2 R^{3} \epsilon \operatorname{a0}^{4}{\left (t \right )} \operatorname{a2}{\left (t \right )} - 3 R^{3} \epsilon \operatorname{a0}^{4}{\left (t \right )} \operatorname{e2}{\left (t \right )} - R^{3} \epsilon \operatorname{a0}^{3}{\left (t \right )} \operatorname{a1}^{2}{\left (t \right )} - 6 R^{3} \epsilon \operatorname{a0}^{3}{\left (t \right )} \operatorname{a1}{\left (t \right )} \operatorname{e1}{\left (t \right )} - 3 R^{3} \epsilon \operatorname{a0}^{3}{\left (t \right )} \operatorname{e1}^{2}{\left (t \right )} + 3 R^{3} \sigma \operatorname{a0}^{2}{\left (t \right )} \operatorname{e2}{\left (t \right )} + 3 R^{3} \sigma \operatorname{a0}{\left (t \right )} \operatorname{e1}^{2}{\left (t \right )} - R^{3} \operatorname{a0}^{3}{\left (t \right )} \frac{d}{d t} \operatorname{a2}{\left (t \right )} - 3 R^{3} \operatorname{a0}^{2}{\left (t \right )} \operatorname{e1}{\left (t \right )} \frac{d}{d t} \operatorname{a1}{\left (t \right )} - 3 R^{3} \operatorname{a0}^{2}{\left (t \right )} \operatorname{e2}{\left (t \right )} \frac{d}{d t} \operatorname{a0}{\left (t \right )} - 3 R^{3} \operatorname{a0}{\left (t \right )} \operatorname{e1}^{2}{\left (t \right )} \frac{d}{d t} \operatorname{a0}{\left (t \right )} - \beta \delta_{P}^{3} k_{1} \operatorname{a2}{\left (t \right )} - \delta_{A} \delta_{P}^{3} k_{1} \operatorname{a2}{\left (t \right )} - 2 \delta_{P}^{3} \epsilon k_{1} \operatorname{a0}{\left (t \right )} \operatorname{a2}{\left (t \right )} - \delta_{P}^{3} \epsilon k_{1} \operatorname{a1}^{2}{\left (t \right )} - \delta_{P}^{3} k_{1} \frac{d}{d t} \operatorname{a2}{\left (t \right )}$$

In [15]:

    




    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-15-3dcb857275fa> in <module>()
----> 1 lm = Matrix([[-(beta + delta_A) - epsilon*A, f1(p)*alpha, R**3*E(t)**2/(delta_P**3*(R**3*E(t)**3/delta_P**3 + k1))],
      2              [8*beta*f2(p), - f1(p)*alpha - delta_M, 0],
      3              [60*beta*(1 - f2(p)), 0, -delta_E]])

NameError: name 'delta_A' is not defined

In [9]:

    

repr(f1(p)/E)


    

    
Out[9]:





'R**3*E(t)**2/(delta_P**3*(R**3*E(t)**3/delta_P**3 + k1))'

In [ ]: