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, naming_scheme=('tau', 'nu', 'kappa'), new_indices_start_at=0)[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.
naming_scheme (tuple) – A tuple of strings giving the names of the variables after reduction. Default is (‘tau’, ‘nu’, ‘kappa’) The 0th element of the tuple is the new name of the independent (time) variable. The 1th is either a string for the pattern for names of independent variables, or a list of strings equal to the number of variables. The 2th element is the pattern for the scaling-invariant constants computed by the algorithm. Since we don’t know how many constants there will be, only a pattern is usable, so the 2th element of the naming_scheme must be a single string.
new_indices_start_at (int) – Where the indices for the new variables should start at. Default is 0. Applies to independent variables if not using a user-supplied list, and always to the reduced constants.

classmethod ODETranslation.from_ode_system(ode_system, naming_scheme=('tau', 'nu', 'kappa'), new_indices_start_at=0)[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, naming_scheme=None, include_aux_vars=True)[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.
naming_scheme – the scheme to use. A tuple of length 3. See

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, include_aux_vars=True)[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

>>> translation.translate_general(system=system, include_aux_vars=False)
dt/dt = 1
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 [HL13].

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, naming_scheme=('t',['n','p'], 'c'))
>>> 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)
{t: t, n: n, p: p, K: 1, d: c0, h: 1, k: c1, r: 1, s: c2}
>>> 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 [HL13]. 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*(1 - t)*exp(t) + t), c1*t*exp(t)*exp(-c1*(1 - t)*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 [HL13]. 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 range(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()
{t: s*t, n: n/d, p: k*p/(d*s), K: K/d, d: 1, h: h*s/k, k: 1, r: r/s, s: 1}

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 [FO99] or [Hub07].

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 [FO99] 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 [FO99] for more details.