Systems of ODEs¶

The ODESystem class is what we use to represent our system of ordinary differential equations.

Creating ODESystems¶

There are a number of different ways of creating an ODESystem.

class desr.ode_system.ODESystem(variables, derivatives, indep_var=None, initial_conditions=None)[source]

A system of differential equations.

The main attributes are variables and derivatives. variables is an ordered tuple of variables, which includes the independent variable. derivatives is an ordered tuple of the same length that contains the derivatives with respect to indep_var.

Parameters:	variables (tuple of sympy.Symbol) – Ordered tuple of variables. derivatives (tuple of sympy.Expression) – Ordered tuple of derivatives. indep_var (sympy.Symbol, optional) – Independent variable we are differentiating with respect to. initial_conditions (tuple of sympy.Symbol) – The initial values of non-constant variables

classmethod ODESystem.from_equations(equations, indep_var=t, initial_conditions=None)[source]

Instantiate from multiple equations.

Parameters:	equations (str, iter of str) – Equations of the form “dx/dt = expr”, optionally seperated by `\n`. indep_var (sympy.Symbol) – The independent variable, usually `t`. initial_conditions (tuple of sympy.Symbol) – The initial values of non-constant variables
Returns:	System of equations.
Return type:	ODESystem

>>> eqns = ['dx/dt = c_0*x*y', 'dy/dt = c_1*(1-x)*(1-y)']
>>> ODESystem.from_equations(eqns)
dt/dt = 1
dx/dt = c_0*x*y
dy/dt = c_1*(-x + 1)*(-y + 1)
dc_0/dt = 0
dc_1/dt = 0
>>> eqns = '\n'.join(['dy/dx = c_0*x*y', 'dz/dx = c_1*(1-y)*z**2'])
>>> ODESystem.from_equations(eqns, indep_var=sympy.Symbol('x'))
dx/dx = 1
dy/dx = c_0*x*y
dz/dx = c_1*z**2*(-y + 1)
dc_0/dx = 0
dc_1/dx = 0

classmethod ODESystem.from_dict(deriv_dict, indep_var=t, initial_conditions=None)[source]

Instantiate from a text of equations.

Parameters:	deriv_dict (dict) – {variable: derivative} mapping. indep_var (sympy.Symbol) – Independent variable, that the derivatives are with respect to. initial_conditions (tuple of sympy.Symbol) – The initial values of non-constant variables
Returns:	System of ODEs.
Return type:	ODESystem

>>> _input = {'x': 'c_0*x*y', 'y': 'c_1*(1-x)*(1-y)'}
>>> _input = {sympy.Symbol(k): sympy.sympify(v) for k, v in _input.iteritems()}
>>> ODESystem.from_dict(_input)
dt/dt = 1
dx/dt = c_0*x*y
dy/dt = c_1*(-x + 1)*(-y + 1)
dc_0/dt = 0
dc_1/dt = 0

>>> _input = {'y':  'c_0*x*y', 'z': 'c_1*(1-y)*z**2'}
>>> _input = {sympy.Symbol(k): sympy.sympify(v) for k, v in _input.iteritems()}
>>> ODESystem.from_dict(_input, indep_var=sympy.Symbol('x'))
dx/dx = 1
dy/dx = c_0*x*y
dz/dx = c_1*z**2*(-y + 1)
dc_0/dx = 0
dc_1/dx = 0

classmethod ODESystem.from_tex(tex)[source]

Given the LaTeX of a system of differential equations, return a ODESystem of it.

Parameters:	tex (str) – LaTeX
Returns:	System of ODEs.
Return type:	ODESystem

>>> eqns = ['\frac{dE}{dt} &= - k_1 E S + k_{-1} C + k_2 C \\',
... '\frac{dS}{dt} &= - k_1 E S + k_{-1} C \\',
... '\frac{dC}{dt} &= k_1 E S - k_{-1} C - k_2 C \\',
... '\frac{dP}{dt} &= k_2 C']
>>> ODESystem.from_tex('\n'.join(eqns))
dt/dt = 1
dC/dt = -C*k_2 - C*k_m1 + E*S*k_1
dE/dt = C*k_2 + C*k_m1 - E*S*k_1
dP/dt = C*k_2
dS/dt = C*k_m1 - E*S*k_1
dk_1/dt = 0
dk_2/dt = 0
dk_m1/dt = 0

Todo

Allow initial conditions to be set from tex.

Finding Scaling Actions¶

ODESystem.power_matrix()[source]

Determine the ‘exponent’ or ‘power’ matrix of the system, denoted by \(K\) in the literature, by gluing together the power matrices of each derivative.

In particular, it concatenates \(K_{\left(\frac{t}{x} \cdot \frac{dx}{dt}\right)}\) for \(x\) in variables, where \(t\) is the independent variable.

>>> eqns = '\n'.join(['ds/dt = -k_1*e_0*s + (k_1*s + k_m1)*c',
... 'dc/dt = k_1*e_0*s - (k_1*s + k_m1 + k_2)*c'])
>>> system = ODESystem.from_equations(eqns)
>>> system.variables
(t, c, s, e_0, k_1, k_2, k_m1)
>>> system.power_matrix()
Matrix([
[1, 1, 1,  1, 1, 1,  1],
[0, 0, 0, -1, 0, 1,  1],
[1, 0, 0,  1, 0, 0, -1],
[0, 0, 0,  1, 1, 0,  0],
[1, 0, 0,  1, 1, 1,  0],
[0, 1, 0,  0, 0, 0,  0],
[0, 0, 1,  0, 0, 0,  1]])

While we get a different answer to the example in the paper, this is just due to choosing our reference exponent in a different way.

Todo

Change the code to agree with the paper.

>>> system.update_initial_conditions({'s': 's_0'})
>>> system.power_matrix()
Matrix([
[1, 1, 1,  1, 1, 1,  1,  0],
[0, 0, 0, -1, 0, 1,  1,  0],
[1, 0, 0,  1, 0, 0, -1,  1],
[0, 0, 0,  1, 1, 0,  0,  0],
[1, 0, 0,  1, 1, 1,  0,  0],
[0, 1, 0,  0, 0, 0,  0,  0],
[0, 0, 1,  0, 0, 0,  1,  0],
[0, 0, 0,  0, 0, 0,  0, -1]])

ODESystem.maximal_scaling_matrix()[source]

Determine the maximal scaling matrix leaving this system invariant.

Returns:	Maximal scaling matrix.
Return type:	sympy.Matrix

>>> eqns = '\n'.join(['ds/dt = -k_1*e_0*s + (k_1*s + k_m1)*c',
... 'dc/dt = k_1*e_0*s - (k_1*s + k_m1 + k_2)*c'])
>>> system = ODESystem.from_equations(eqns)
>>> system.maximal_scaling_matrix()
Matrix([
[1, 0, 0, 0, -1, -1, -1],
[0, 1, 1, 1, -1,  0,  0]])

Output Functions¶

There are a number of useful ways to output the system.

ODESystem.derivative_dict

expr mapping, filtering out the None’s in expr.

Returns:	Keys are non-constant variables, value is the derivative with respect to the independent variable.
Return type:	dict
Type:	Return a variable

ODESystem.to_tex()[source]

Returns:	TeX representation.
Return type:	str

>>> eqns = ['dC/dt = -C*k_2 - C*k_m1 + E*S*k_1',
... 'dE/dt = C*k_2 + C*k_m1 - E*S*k_1',
... 'dP/dt = C*k_2',
... 'dS/dt = C*k_m1 - E*S*k_1']
>>> system = ODESystem.from_equations('\n'.join(eqns))
>>> print system.to_tex()
\frac{dt}{dt} &= 1 \\
\frac{dC}{dt} &= - C k_{2} - C k_{-1} + E S k_{1} \\
\frac{dE}{dt} &= C k_{2} + C k_{-1} - E S k_{1} \\
\frac{dP}{dt} &= C k_{2} \\
\frac{dS}{dt} &= C k_{-1} - E S k_{1} \\
\frac{dk_{1}}{dt} &= 0 \\
\frac{dk_{2}}{dt} &= 0 \\
\frac{dk_{-1}}{dt} &= 0

>>> system.update_initial_conditions({'C': 'C_0'})
>>> print system.to_tex()
\frac{dt}{dt} &= 1 \\
\frac{dC}{dt} &= - C k_{2} - C k_{-1} + E S k_{1} \\
\frac{dE}{dt} &= C k_{2} + C k_{-1} - E S k_{1} \\
\frac{dP}{dt} &= C k_{2} \\
\frac{dS}{dt} &= C k_{-1} - E S k_{1} \\
\frac{dk_{1}}{dt} &= 0 \\
\frac{dk_{2}}{dt} &= 0 \\
\frac{dk_{-1}}{dt} &= 0 \\
\frac{dC_{0}}{dt} &= 0 \\
C\left(0\right) &= C_{0}

>>> system.add_constraints('K_m', '(k_m1 + k_2) / k_1')
>>> print system.to_tex()
\frac{dt}{dt} &= 1 \\
\frac{dC}{dt} &= - C k_{2} - C k_{-1} + E S k_{1} \\
\frac{dE}{dt} &= C k_{2} + C k_{-1} - E S k_{1} \\
\frac{dP}{dt} &= C k_{2} \\
\frac{dS}{dt} &= C k_{-1} - E S k_{1} \\
\frac{dk_{1}}{dt} &= 0 \\
\frac{dk_{2}}{dt} &= 0 \\
\frac{dk_{-1}}{dt} &= 0 \\
\frac{dC_{0}}{dt} &= 0 \\
\frac{dK_{m}}{dt} &= 0 \\
C\left(0\right) &= C_{0} \\
K_{m} &= \frac{1}{k_{1}} \left(k_{2} + k_{-1}\right)