ODETranslation and Reduction of ODESystems¶

The ODETranslation is the class used to represent the scaling symmetries of a system.

Creating ODE Translations¶

class desr.ode_translation.ODETranslation(scaling_matrix, variables_domain=None, hermite_multiplier=None)[source]

An object used for translating between systems of ODEs according to a scaling matrix. The key data are scaling_matrix and herm_mult, which together contain all the information needed (along with a variable order) to reduce an ODESystem.

Parameters:

scaling_matrix (sympy.Matrix) – Matrix that defines the torus action on the system.
variables_domain (iter of sympy.Symbol, optional) – An ordering of the variables we expect to act upon. If this is not given, we will act on a system according to the position of the variables in variables, as long as there is the correct number of variables.
hermite_multiplier (sympy.Matrix, optional) – User-defined Hermite multiplier, that puts scaling_matrix into column Hermite normal form. If not given, the normal Hermite multiplier will be calculated.

classmethod ODETranslation.from_ode_system(ode_system)[source]

Create a ODETranslation given an ode_system.ODESystem instance, by taking the maximal scaling matrix.

Parameters:	ode_system (ODESystem) –
Return type:	ODETranslation

Reduction Methods¶

ODETranslation.translate(system)[source]

Translate an ode_system.ODESystem into a reduced system.

First, try the simplest parameter reduction method, then the dependent variable translation (where the scaling action ignores the independent variable) and finally the general reduction scheme.

Parameters:	system (ODESystem) – System to reduce.
Return type:	ODESystem

ODETranslation.translate_parameter(system)[source]: Translate according to parameter scheme

ODETranslation.translate_dep_var(system)[source]

Given a system of ODEs, translate the system into a simplified version. Assume we are only working on dependent variables, not the independent variable.

Parameters:	system (ODESystem) – System to reduce.
Return type:	ODESystem

>>> equations = 'dz1/dt = z1*(1+z1*z2);dz2/dt = z2*(1/t - z1*z2)'.split(';')
>>> system = ODESystem.from_equations(equations)
>>> translation = ODETranslation.from_ode_system(system)
>>> translation.translate_dep_var(system=system)
dt/dt = 1
dx0/dt = x0*(y0 - 1/t)
dy0/dt = y0*(1 + 1/t)

ODETranslation.translate_general(system)[source]

The most general reduction scheme. If there are \(n\) variables (including the independent variable) then there will be a system of \(n-r+1\) invariants and \(r\) auxiliary variables.

Parameters:	system (ODESystem) – System to reduce.
Return type:	ODESystem

>>> equations = 'dn/dt = n*( r*(1 - n/K) - k*p/(n+d) );dp/dt = s*p*(1 - h*p / n)'.split(';')
>>> system = ODESystem.from_equations(equations)
>>> translation = ODETranslation.from_ode_system(system)
>>> translation.translate_general(system=system)
dt/dt = 1
dx0/dt = 0
dx1/dt = 0
dx2/dt = 0
dy0/dt = y0/t
dy1/dt = y1*(-y0*y1*y5/y3 - y0*y2/(y1 + 1) + y0*y5)/t
dy2/dt = y2*(y0 - y0*y2*y4/y1)/t
dy3/dt = 0
dy4/dt = 0
dy5/dt = 0

Reverse Translation¶

Reverse translation is the process of taking solutions of the reduced system and recovering solutions of the original system.

ODETranslation.reverse_translate(variables)[source]

Given the solutions of a reduced system, reverse translate them into solutions of the original system.

Note

This doesn’t reverse translate the parameter reduction scheme.

It only guesses between reverse_translate_general() and reverse_translate_dep_var()

Parameters:	variables (iter of sympy.Expression) – The solution auxiliary variables \(x(t)\) and solution invariants \(y(t)\) of the reduced system. Variables should be ordered auxiliary variables followed by invariants. indep_var_index (int) – The location of the independent variable.
Return type:	tuple

ODETranslation.reverse_translate_parameter(variables)[source]

Given the solutions of a reduced system, reverse translate them into solutions of the original system.

Parameters:	variables (iter of sympy.Expression) – \(r\) constants (auxiliary variables) followed by the independent, dependent and invariant constants.

Example 7.5 (Prey-predator model) from [Hubert2013c].

Use the matrices from the paper, which differ to ours as different conventions are used.

\[ \begin{align}\begin{aligned}\frac{dn}{dt} = n \left( r \left(1 - \frac{n}{K} \right) - \frac{k p}{n + d} \right)\\\frac{dp}{dt} = s p \left(1 - \frac{h p}{n} \right)\end{aligned}\end{align} \]

>>> equations = ['dn/dt = n*( r*(1 - n/K) - k*p/(n+d) )','dp/dt = s*p*(1 - h*p / n)']
>>> system = ODESystem.from_equations(equations)
>>> system.variables
(t, n, p, K, d, h, k, r, s)
>>> maximal_scaling_matrix = system.maximal_scaling_matrix()
>>> # Hand craft a Hermite multiplier given the invariants from the paper.
>>> # This is done by placing the auxiliaries first, then rearranging the invariants and reading them off
>>> # individually from the Hermite multiplier given in the paper
>>> herm_mult = sympy.Matrix([[ 0, 0,  0,  1,  0,  0,  0,  0,  0],  # t
...                           [ 0, 0,  0,  0,  1,  0,  0,  0,  0],  # n
...                           [ 0, 0,  0,  0,  0,  1,  0,  0,  0],  # p
...                           [ 0, 1,  1,  0, -1, -1, -1,  0,  0],  # K
...                           [ 0, 0,  0,  0,  0,  0,  1,  0,  0],  # d
...                           [ 0, 0, -1,  0,  0,  1,  0, -1,  0],  # h
...                           [ 0, 0,  0,  0,  0,  0,  0,  1,  0],  # k
...                           [-1, 0,  0,  1,  0,  0,  0, -1, -1],  # r
...                           [ 0, 0,  0,  0,  0,  0,  0,  0,  1]])  # s
>>> translation = ODETranslation(maximal_scaling_matrix,
...                              variables_domain=system.variables,
...                              hermite_multiplier=herm_mult)
>>> translation.translate_parameter(system)
dt/dt = 1
dn/dt = -c1*n*p/(c0 + n) - n**2 + n
dp/dt = c2*p - c2*p**2/n
dc0/dt = 0
dc1/dt = 0
dc2/dt = 0
>>> translation.translate_parameter_substitutions(system)
{k: c1, n: n, r: 1, d: c0, K: 1, h: 1, s: c2, p: p, t: t}
>>> soln_reduced = sympy.var('x1, x2, x3, t, n, p, c0, c1, c2')
>>> translation.reverse_translate_parameter(soln_reduced)
Matrix([[t*x1, n*x2, p*x3, x2, c0*x2, x2/x3, c1*x2/(x1*x3), 1/x1, c2/x1]])

ODETranslation.reverse_translate_dep_var(variables, indep_var_index)[source]

Given the solutions of a (dep_var) reduced system, reverse translate them into solutions of the original system.

Parameters:	variables (iter of sympy.Expression) – The solution auxiliary variables \(x(t)\) and solution invariants \(y(t)\) of the reduced system. Variables should be ordered auxiliary variables followed by invariants. indep_var_index (int) – The location of the independent variable.
Return type:	tuple

Example 6.4 from [Hubert2013c]. We will just take the matrices from the paper, rather than use our own (which are constructed using different conventions).

\[ \begin{align}\begin{aligned}\frac{dz_1}{dt} = z_1 \left( 1+ z_1 z_2 \right)\\\frac{dz_2}{dt} = z_2 \left(\frac{1}{t} - z_1 z_2 \right)\end{aligned}\end{align} \]

>>> equations = 'dz1/dt = z1*(1+z1*z2);dz2/dt = z2*(1/t - z1*z2)'.split(';')
>>> system = ODESystem.from_equations(equations)
>>> var_order = ['t', 'z1', 'z2']
>>> system.reorder_variables(var_order)
>>> scaling_matrix = sympy.Matrix([[0, 1, -1]])
>>> hermite_multiplier = sympy.Matrix([[0, 1, 0],
...                                    [1, 0, 1],
...                                    [0, 0, 1]])
>>> translation = ODETranslation(scaling_matrix=scaling_matrix,
...                              hermite_multiplier=hermite_multiplier)
>>> reduced = translation.translate_dep_var(system)
>>> # Check reduction
>>> answer = ODESystem.from_equations('dx0/dt = x0*(y0 + 1);dy0/dt = y0*(1 + 1/t)'.split(';'))
>>> reduced == answer
True
>>> reduced
dt/dt = 1
dx0/dt = x0*(y0 + 1)
dy0/dt = y0*(1 + 1/t)
>>> t_var, c1_var, c2_var = sympy.var('t c1 c2')
>>> reduced_soln = (c2_var*sympy.exp(t_var+c1_var*(1-t_var)*sympy.exp(t_var)),
...                 c1_var * t_var * sympy.exp(t_var))
>>> original_soln = translation.reverse_translate_dep_var(reduced_soln, system.indep_var_index)
>>> original_soln == (reduced_soln[0], reduced_soln[1] / reduced_soln[0])
True
>>> original_soln
(c2*exp(c1*(-t + 1)*exp(t) + t), c1*t*exp(t)*exp(-c1*(-t + 1)*exp(t) - t)/c2)

ODETranslation.reverse_translate_general(variables, system_indep_var_index=0)[source]

Given an iterable of variables, or exprs, reverse translate into the original variables. Here we expect t as the first variable, since we need to divide by it and substitute

Given the solutions of a (general) reduced system, reverse translate them into solutions of the original system.

Parameters:	variables (iter of sympy.Expression) – The independent variable \(t\), the solution auxiliary variables \(x(t)\) and the solution invariants \(y(t)\) of the reduced system. Variables should be ordered: independent variable, auxiliary variables then invariants. indep_var_index (int) – The location of the independent variable.
Return type:	tuple

Example 6.6 from [Hubert2013c]. We will just take the matrices from the paper, rather than use our own (which are constructed using different conventions).

\[ \begin{align}\begin{aligned}\frac{dz_1}{dt} = \frac{z_1 \left(z_1^5 z_2 - 2 \right)}{3t}\\\frac{dz_2}{dt} = \frac{z_2 \left(10 - 2 z_1^5 z_2 + \frac{3 z_1^2 z_2}{t} \right)}{3t}\end{aligned}\end{align} \]

>>> equations = ['dz1/dt = z1*(z1**5*z2 - 2)/(3*t)',
...              'dz2/dt = z2*(10 - 2*z1**5*z2 + 3*z1**2*z2/t )/(3*t)']
>>> system = ODESystem.from_equations(equations)
>>> var_order = ['t', 'z1', 'z2']
>>> system.reorder_variables(var_order)
>>> scaling_matrix = sympy.Matrix([[3, -1, 5]])
>>> hermite_multiplier = sympy.Matrix([[1, 1, -1],
...                                    [2, 3, 2],
...                                    [0, 0, 1]])
>>> translation = ODETranslation(scaling_matrix=scaling_matrix,
...                              hermite_multiplier=hermite_multiplier)
>>> reduced = translation.translate_general(system)
>>> # Quickly check the reduction
>>> answer = ODESystem.from_equations(['dx0/dt = x0*(2*y0*y1/3 - 1/3)/t',
...                                    'dy0/dt = y0*(y0*y1 - 1)/t',
...                                    'dy1/dt = y1*(y1 + 1)/t'])
>>> reduced == answer
True
>>> reduced
dt/dt = 1
dx0/dt = x0*(2*y0*y1/3 - 1/3)/t
dy0/dt = y0*(y0*y1 - 1)/t
dy1/dt = y1*(y1 + 1)/t

Check reverse translation:

>>> reduced_soln = (sympy.var('t'),
...                 sympy.sympify('c3/(t**(1/3)*(ln(t-c1)-ln(t)+c2)**(2/3))'),  # x
...                 sympy.sympify('c1/(t*(ln(t-c1)-ln(t)+c2))'),  # y1
...                 sympy.sympify('t/(c1 - t)'))  # y2
>>> original_soln = translation.reverse_translate_general(reduced_soln, system.indep_var_index)
>>> original_soln_answer = [#reduced_soln[1] ** 3 / (reduced_soln[2] ** 2)  # x^3 / y1^2
...                         reduced_soln[2] / reduced_soln[1],  # y1 / x
...                         reduced_soln[1] ** 5 * reduced_soln[3] / reduced_soln[2] ** 4]  # x^5 y2 / y1^4
>>> original_soln_answer = tuple([soln.subs({sympy.var('t'): sympy.sympify('t / (c3**3 / c1**2)')})
...                               for soln in original_soln_answer])
>>> original_soln == original_soln_answer
True
>>> original_soln[0]
c1/(c3*(c1**2*t/c3**3)**(2/3)*(c2 - log(c1**2*t/c3**3) + log(c1**2*t/c3**3 - c1))**(1/3))
>>> original_soln[1]
c3**5*(c1**2*t/c3**3)**(10/3)*(c2 - log(c1**2*t/c3**3) + log(c1**2*t/c3**3 - c1))**(2/3)/(c1**4*(-c1**2*t/c3**3 + c1))

Extending from a set of Invariants¶

ODETranslation.extend_from_invariants(invariant_choice)[source]

Extend a given set of invariants, expressed as a matrix of exponents, to find a Hermite multiplier that will rewrite the system in terms of invariant_choice.

Parameters:	invariant_choice (sympy.Matrix) – The \(n \times k\) matrix representing the invariants, where \(n\) is the number of variables and \(k\) is the number of given invariants.
Returns:	An ODETranslation representing the rewrite rules in terms of the given invariants.
Return type:	ODETranslation

>>> variables = sympy.symbols(' '.join(['y{}'.format(i) for i in xrange(6)]))
>>> ode_translation = ODETranslation(sympy.Matrix([[1, 0, 3, 0, 2, 2],
...                                                [0, 2, 0, 1, 0, 1],
...                                                [2, 0, 0, 3, 0, 0]]))
>>> ode_translation.invariants(variables=variables)
Matrix([[y0**3*y2*y5**2/(y3**2*y4**5), y1*y4**2/y5**2, y2**2/y4**3]])

Now we can express two new invariants, which we think are more interesting, as a matrix. We pick the product of the first two invariants, and the product of the last two invariants: y0**3 * y1 * y2/(y3**2 * y4**3) and y1 * y2**2/(y4 * y5**2)

>>> new_inv = sympy.Matrix([[3, 1, 1, -2, -3, 0],
...                         [0, 1, 2, 0, -1, -2]]).T

>>> new_translation = ode_translation.extend_from_invariants(new_inv)
>>> new_translation.invariants(variables=variables)
Matrix([[y0**3*y1*y2/(y3**2*y4**3), y1*y2**2/(y4*y5**2), y1*y4**2/y5**2]])

Useful Attributes¶

ODETranslation.invariants(variables=None)[source]

The invariants of the system, as determined by herm_mult_n.

Returns:	The invariants of the system, in terms of `variables_domain` if not `None`. Otherwise, use `variables` and failing that use \(y_i\).
Return type:	tuple

>>> equations = ['dn/dt = n*( r*(1 - n/K) - k*p/(n+d) )','dp/dt = s*p*(1 - h*p / n)']
>>> system = ODESystem.from_equations(equations)
>>> translation = ODETranslation.from_ode_system(system)
>>> translation.invariants()
Matrix([[s*t, n/d, k*p/(d*s), K/d, h*s/k, r/s]])

ODETranslation.auxiliaries(variables=None)[source]

The auxiliary variables of the system, as determined by herm_mult_i.

Returns:	The auxiliary variables of the system, in terms of `variables_domain` if not `None`. Otherwise, use `variables` and failing that use \(x_i\).
Return type:	tuple

>>> equations = ['dn/dt = n*( r*(1 - n/K) - k*p/(n+d) )','dp/dt = s*p*(1 - h*p / n)']
>>> system = ODESystem.from_equations(equations)
>>> translation = ODETranslation.from_ode_system(system)
>>> translation.auxiliaries()
Matrix([[1/s, d, d*s/k]])

ODETranslation.rewrite_rules(variables=None)[source]

Given a set of variables, print the rewrite rules. These are the rules that allow you to write any rational invariant in terms of the invariants.

Parameters:	variables (iter of sympy.Symbol) – The names of the variables appearing in the rational invariant.
Returns:	A `dict` whose keys are the given variables and whose values are the values to be substituted.
Return type:	dict

>>> equations = ['dn/dt = n*( r*(1 - n/K) - k*p/(n+d) )','dp/dt = s*p*(1 - h*p / n)']
>>> system = ODESystem.from_equations(equations)
>>> translation = ODETranslation.from_ode_system(system)
>>> translation.rewrite_rules()
{k: 1, n: n/d, r: r/s, d: 1, K: K/d, h: h*s/k, s: 1, p: k*p/(d*s), t: s*t}

For example, \(r t\) is an invariant. Substituting in the above mapping gives us a way to write it in terms of our generating set of invariants:

\[r t \mapsto \left( \frac{r}{s} \right) \left( s t \right) = r t\]

Output Functions¶

ODETranslation.to_tex()[source]

Returns:	The scaling matrix \(A\), the Hermite multiplier \(V\) and \(W = V^{-1}\), in beautiful LaTeX.
Return type:	str

>>> print ODETranslation(sympy.Matrix(range(12)).reshape(3, 4)).to_tex()
A=
0 & 1 & 2 & 3 \\
4 & 5 & 6 & 7 \\
8 & 9 & 10 & 11 \\
V=
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
-1 & 3 & -3 & -2 \\
1 & -2 & 2 & 1 \\
W=
0 & 1 & 2 & 3 \\
1 & 1 & 1 & 1 \\
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\

Advanced Methods¶

These methods will be familiar to those who use Lie groups to analyse more general symmetries of differential equations. For more information, see [Fels1999] or [Hubert2007].

ODETranslation.moving_frame(variables=None)[source]

Given a finite-dimensional Lie group \(G\) acting on a manifold \(M\), a moving frame is defined as a \(G\)-equivariant map \(\rho : M \rightarrow G\).

We can construct a moving frame given a rational section, by sending a point \(x \in M\) to the Lie group element \(\rho ( x )\) such that \(\rho ( x ) \cdot x\) lies on \(K\).

>>> equations = 'dz1/dt = z1*(1+z1*z2);dz2/dt = z2*(1/t - z1*z2)'.split(';')
>>> system = ODESystem.from_equations(equations)
>>> translation = ODETranslation.from_ode_system(system)
>>> moving_frame = translation.moving_frame()
>>> moving_frame
(z2,)

Now we can check that the moving frame moves a point to the cross-section.

>>> translation.rational_section()
Matrix([[1 - 1/z2]])
>>> scale_action(moving_frame, translation.scaling_matrix)  # \rho(x)
Matrix([[1, z2, 1/z2]])
>>> # :math:`\rho(x).x` should satisfy the equations of the rational section.
>>> scale_action(moving_frame, translation.scaling_matrix).multiply_elementwise(sympy.Matrix(system.variables).T)
Matrix([[t, z1*z2, 1]])

Another example:

>>> equations = ['dn/dt = n*( r*(1 - n/K) - k*p/(n+d) )','dp/dt = s*p*(1 - h*p / n)']
>>> system = ODESystem.from_equations(equations)
>>> translation = ODETranslation.from_ode_system(system)
>>> moving_frame = translation.moving_frame()
>>> moving_frame
(s, 1/d, k/(d*s))

Now we can check that the moving frame moves a point to the cross-section.

>>> translation.rational_section()
Matrix([[1 - 1/s, d - 1, d*s - 1/k]])
>>> scale_action(moving_frame, translation.scaling_matrix)  # \rho(x)
Matrix([[s, 1/d, k/(d*s), 1/d, 1/d, s/k, 1/k, 1/s, 1/s]])
>>> # :math:`\rho(x).x` should satisfy the equations of the rational section.
>>> scale_action(moving_frame, translation.scaling_matrix).multiply_elementwise(sympy.Matrix(system.variables).T)
Matrix([[s*t, n/d, k*p/(d*s), K/d, 1, h*s/k, 1, r/s, 1]])

See [Fels1999] for more details.

ODETranslation.rational_section(variables=None)[source]

Given a finite-dimensional Lie group \(G\) acting on a manifold \(M\), a cross-section \(K \subset M\) is a subset that intersects each orbit of \(M\) in exactly one place.

Parameters:	variables (iter of sympy.Symbol, optional) – Variables of the system.
Returns:	A set of equations defining a variety that is a cross-section.
Return type:	sympy.Matrix

>>> equations = 'dz1/dt = z1*(1+z1*z2);dz2/dt = z2*(1/t - z1*z2)'.split(';')
>>> system = ODESystem.from_equations(equations)
>>> translation = ODETranslation.from_ode_system(system)
>>> translation.rational_section()
Matrix([[1 - 1/z2]])

>>> equations = ['dn/dt = n*( r*(1 - n/K) - k*p/(n+d) )','dp/dt = s*p*(1 - h*p / n)']
>>> system = ODESystem.from_equations(equations)
>>> translation = ODETranslation.from_ode_system(system)
>>> translation.rational_section()
Matrix([[1 - 1/s, d - 1, d*s - 1/k]])

See [Fels1999] for more details.