Available Transformations

Carleman linearization of polynomial differential equations in SageMath.

Reduction methods

quadratic_reduction() Reduce from any order in standard form to quadratic order
transfer_matrices() Compute the higher order transfer matrices \(A^{i}_{i+j-1}\)

Linearization

linearize() Compute Carleman linearization and export to a MAT file
truncated_matrix() Finite order \(N\) Carleman linearization

Kronecker power and linear algebra

get_index_from_key() Return first occurrence of given key over a Kronecker power
get_key_from_index() Return multi-index of Kronecker power given an index and the order
kron_prod() Compute the Kronecker product of x and y
kron_power() Receives a \(n\times 1\) vector and computes its Kronecker power \(x^{[i]}\)
log_norm() Compute the logarithmic norm of a matrix

Error computation

characteristics() Information about the norms of the matrices in \(F\)
error_function() Compute the error function of a truncated ODE
plot_error_function() Plot the estimated error of the linearized ODE as a function of time

AUTHOR:

  • Marcelo Forets (Dec 2016 at VERIMAG - France)

MF acknowledges the hospitality at Max-Planck Institute for Software Systems, Saarbrucken, Germany, where part of this package was written (Apr 2016).

carlin.transformation.characteristics(F, n, k, ord=inf)

Information about the norms of the matrices in \(F\).

INPUT:

  • F – list of matrices in some NumPy sparse format, for which the toarray method is available
  • n – dimension on state-space
  • k – order of the system
  • ord – order of the \(p\)-th norm, for \(1 \leq p < \infty\), and p='inf' for \(p=\infty\)

OUTPUT:

Dictionary c containing norm_Fi_inf, log_norm_F1_inf and beta0_const.

carlin.transformation.error_function(model, N, x0)

Compute the error function of a truncated ODE.

INPUT:

  • model – Polynomial ODE or string containing the model in text format
  • N – integer; truncation order
  • x0 – list; initial point

OUTPUT:

  • Ts – convergence time computed from the reduced quadratic system
  • error – function of \(t\), the estimated truncation error in the supremum norm

EXAMPLES:

sage: from carlin.transformation import error_function
sage: from carlin.library import quadratic_scalar as P
sage: Ts, error = error_function(P(0.5, 2), 2, [0, 0.5])
sage: Ts
0.8109...
sage: error
0.5*(2.0*e^(0.5*t) - 2.0)^2*e^(0.5*t)/(-2.0*e^(0.5*t) + 3.0)
carlin.transformation.linearize(model, target_filename, N, x0, **kwargs)

Compute Carleman linearization and export to a MAT file.

INPUT:

  • mode – model as a PolynomialODE or a string containing the model in text format
  • target_filename – string with the name of the output file in MAT format
  • N – truncation order
  • x0 – initial point, can be either a list or a polyhedron; see the code for further details

NOTES:

This function is self-contained; it transforms to canonical quadratic form, then computes Carleman linearization together with the error estimates and exports the resulting matrix \(A_N\) and characteristics to a MAT file.

carlin.transformation.plot_error_function(model_filename, N, x0, Tfrac=0.8)

Plot the estimated error of the linearized ODE as a function of time.

INPUT:

  • model_filename – string containing the model in text format
  • N – truncation order
  • x0 – initial point, a list
  • Tfrac – (optional, default: \(0.8\)): fraction of the convergence radius, to specify the plotting range in the time axis

NOTE:

This function calls error_function for the error computations.

carlin.transformation.quadratic_reduction(F, n, k)

Reduce a \(k\)-th order system of polynomial ODE’s into a quadratic one.

INPUT:

  • F – list of matrices defining the system of polynomial ODE’s
  • n – integer, system’s dimension
  • k – integer, order of the polynomial ODE

OUTPUT:

The list [F_tilde, nquad, kquad] corresponding to the reduced list of \(F_j\)‘s, dimension and order respectively.

EXAMPLES:

Consider the following two-dimensional system:

sage: from carlin.library import quadratic_scalar
sage: P = quadratic_scalar(1, -1); P.funcs()
[-x0^2 + x0]
sage: from carlin.io import get_Fj_from_model
sage: (F, n, k) = get_Fj_from_model(P.funcs(), P.dim(), P.degree())
sage: (n, k)
(1, 2)
sage: [matrix(Fi.toarray()) for Fi in F]
[[1.0], [-1.0]]

Since it is already quadratic, the reduction does nothing:

sage: from carlin.transformation import quadratic_reduction
sage: (Fred, nred, kred) = quadratic_reduction(F, n, k)
sage: nred, kred
(1, 2)
sage: [matrix(Fi.toarray()) for Fi in F]
[[1.0], [-1.0]]

Now consider the more interesting case of a cubic system:

sage: from carlin.library import cubic_scalar
sage: P = cubic_scalar(1, -1); P.funcs()
[-x0^3 + x0]
sage: from carlin.io import get_Fj_from_model
sage: (F, n, k) = get_Fj_from_model(P.funcs(), P.dim(), P.degree())
sage: (n, k)
(1, 3)
sage: [matrix(Fi.toarray()) for Fi in F]
[[1.0], [0.0], [-1.0]]

Introducing the auxiliary variables \(\tilde{x}_1 := x\) and \(\tilde{x}_2:=x^2\), the corresponding quadratic system in \(\tilde{x} := (\tilde{x}_1, \tilde{x}_2)\) is:

sage: (Fred, nred, kred) = quadratic_reduction(F, n, k)
sage: nred, kred
(2, 2)
sage: matrix(Fred[0].toarray())
[1.0  0.0]
[0.0  2.0]
sage: matrix(Fred[1].toarray())
[ 0.0 -1.0  0.0  0.0]
[ 0.0  0.0  0.0 -2.0]
carlin.transformation.transfer_matrices(N, F, n, k)

Higher order transfer matrices \(A^{i}_{i+j-1}\).

INPUT:

  • N – order of truncation
  • F – list, sequence of matrices \(F_j\)
  • n – the dimension of the state-space
  • k – the order of the polynomial vector field. It is equal to len(F)

OUTPUT:

  • A – the transfer matrices \(A^{i}_{i+j-1}\) that correspond to \(i = 1, \ldots, N\).

    It is given as a list of lists. Each inner list has dimension \(k\).

carlin.transformation.truncated_matrix(N, *args, **kwargs)

Finite order Carleman linearization.

INPUT:

  • N – order of truncation
  • input_format – valid options are:
    • 'model_filename' – (default); file in text format
    • 'transfer_matrices' – for \((A, n, k)\), which should be given separately
    • 'Fj_matrices' – for \((F, n, k)\), which should be given separately

OUTPUT:

The transfer matrices \(A^{i}_{i+j-1}\) that correspond to \(i = 1, \ldots, N\). It is given as a list of lists, and each inner list has dimension \(k\).

EXAMPLES:

Conisder the polynomial ODE:

sage: from carlin.polynomial_ode import PolynomialODE
sage: x = polygens(QQ, ["x0", "x1"])
sage: f = [x[0]^3*x[1], -2*x[0]+2*x[1]^2]
sage: P = PolynomialODE(f, 2, 4)

Compute the Carleman matrix arising from linearization at order \(N=2\):

sage: from carlin.transformation import get_Fj_from_model, truncated_matrix
sage: Fj = get_Fj_from_model(P.funcs(), P.dim(), P.degree())
sage: matrix(truncated_matrix(2, *Fj, input_format="Fj_matrices").toarray())
[ 0.0  0.0  0.0  0.0  0.0  0.0]
[-2.0  0.0  0.0  0.0  0.0  2.0]
[ 0.0  0.0  0.0  0.0  0.0  0.0]
[ 0.0  0.0 -2.0  0.0  0.0  0.0]
[ 0.0  0.0 -2.0  0.0  0.0  0.0]
[ 0.0  0.0  0.0 -2.0 -2.0  0.0]

Try a higher truncation order:

sage: matrix(truncated_matrix(4, *Fj, input_format="Fj_matrices").toarray())
30 x 30 dense matrix over Real Double Field (use the '.str()' method to see the entries)

Handling Input/Output

Carleman linearization input/output functions.

Load and export polynomial ODEs

export_model_to_mat() Export model to a MAT file as the sequence of sparse \(F_j\) matrices
load_model() Read an ODE system from a text file

Transformation functions

get_Fj_from_model() Transform a model into standard form as a sum of Kronecker products

Solve polynomial ODEs

solve_ode_exp() Solve Carleman linearized ODE.
plot_truncated() Solve and return the graphics in phase space.

AUTHOR:

  • Marcelo Forets (Dec 2016 at VERIMAG - UGA)
carlin.io.export_model_to_mat(model_filename, F=None, n=None, k=None, **kwargs)

Export ODE model to a Matlab .mat format.

INPUT:

The model can be given either as a model in text file, or as the tuple \((F, n, k)\).

  • model_filename – string with the model filename. If \((F, n, k)\) is not provided, then such model should be reachable from the current path
  • F – list of sparse matrices \(F_1, \ldots, F_k\)
  • n – dimension of state-space
  • k – order of the system

OUTPUT:

List containing the solution of the 1st order ODE \(x'(t) = A_N x(t)\), with initial condition \(x(0) = x_0\). The output filename is model_filename, replacing the .sage extension with the .mat extension.

carlin.io.get_Fj_from_model(model_filename=None, f=None, n=None, k=None)

Transform an input model of a polynomial vector field into standard form as a sum of Kronecker products.

The model can be given either as an external file (model_filename), or as the data \(f\), \(n\) and \(k\).

INPUT:

  • model_filename – string containing the filename

OUTPUT:

  • F – list of sparse matrices \(F_1, \ldots, F_k\). These are formatted in dictionary-of-keys (dok) form.
  • n – dimension of the state-space of the system
  • k – degree of the system

EXAMPLES:

Let’s create a two-dimensional system of second order:

sage: x = polygens(QQ, ['x0', 'x1'])
sage: f = [-2*x[0] + x[1] + x[0]**2, x[0] + 3/2*x[0]*x[1]]
sage: from carlin.polynomial_ode import PolynomialODE
sage: T = PolynomialODE(f, 2, 2); T
A Polynomial ODE in n = 2 variables

There are two \(F_j\) matrices in sparse representation, which can be visualizeda with the toarray method:

 sage: from carlin.io import get_Fj_from_model
 sage: F, _, _ = get_Fj_from_model(T.funcs(), T.dim(), T.degree())
 sage: F
 [<2x2 sparse matrix of type '<type 'numpy.float64'>'
         with 3 stored elements in Dictionary Of Keys format>,
  <2x4 sparse matrix of type '<type 'numpy.float64'>'
         with 2 stored elements in Dictionary Of Keys format>]
 sage: F[0].toarray()
 array([[-2.,  1.],
        [ 1.,  0.]])
sage: F[1].toarray()
array([[ 1. ,  0. ,  0. ,  0. ],
       [ 0. ,  1.5,  0. ,  0. ]])

TESTS:

Check if a polynomial with a non-numeric coefficient is allowed:

sage: from carlin.io import get_Fj_from_model sage: x = polygens(QQ, [“x0”]); mu = SR.var(“mu”) sage: from carlin.polynomial_ode import PolynomialODE sage: T = PolynomialODE([mu*x[0]^2 - x[0]], 1, 2) sage: get_Fj_from_model(T.funcs(), T.dim(), T.degree()) Traceback (most recent call last): ... NotImplementedError: the coefficients of the polynomial should be numeric

carlin.io.load_model(model_filename)

Read an input ODE system.

INPUT:

  • model_filename – string with the model filename

OUTPUT:

  • f – list of multivariate polynomials which describes the system of ODEs, each component in the polynomial ring \(\mathbb{Q}[x_1,\ldots,x_n]\)
  • n – integer, dimension of f
  • k – integer, degree of f
carlin.io.plot_truncated(model, N, x0, tini, T, NPOINTS, xcoord=0, ycoord=1, **kwargs)

Solve and return graphics in phase space of a given model.

INPUT:

  • model – PolynomialODE, defining the tuple \((f, n, k)\)
  • N – integer; truncation order
  • x0 – vector; initial condition
  • tini – initial time of simulation
  • T – final time of simulation
  • NPOINTS – number of points sampled
  • xcoord – (default: \(0\)), x-coordinate in plot
  • ycoord – (default: \(1\)), y coordinate in plot

NOTES:

By default, returns a plot in the plane \((x_1, x_2)\). All other keyword arguments passes are sent to the \(list_plot\) command (use to set line color, style, etc.)

EXAMPLES:

sage: from carlin.library import vanderpol
sage: from carlin.io import plot_truncated
sage: G = plot_truncated(vanderpol(1, 1), 2, [0.1, 0], 0, 5, 100)
sage: G.show(gridlines=True, axes_labels = ['$x_1$', '$x_2$']) # not tested
Graphics object consisting in 1 graphics primitive

All other keyword arguments are passed to the list_plot function. For example, specify color and maximum and minimum values for the axes:

sage: G = plot_truncated(vanderpol(1, 1), 2, [0.1, 0], 0, 5, 100, color='green', xmin=-1, xmax=1, ymin=-1, ymax=1)
sage: G.show(gridlines=True, axes_labels = ['$x_1$', '$x_2$']) # not tested
Graphics object consisting in 1 graphics primitive
carlin.io.solve_ode_exp(AN, x0, N, tini=0, T=1, NPOINTS=100)

Solve 1st order ODE with initial condition in Kronecker power form.

INPUT:

  • AN – matrix, it can be Sage dense or NumPy sparse in COO format
  • x0 – vector, initial point
  • N – integer, order of the truncation
  • tini – (optional, default: 0) initial time
  • T – (optional, default: 1) final time
  • NPOINTS – (optional, default: 100) number of points computed

OUTPUT:

List containing the solution of the 1st order ODE \(x'(t) = A_N x(t)\), with initial condition \(x(0) = x_0\). The solution is computed the matrix exponential \(e^{A_N t_i}\) directly.

NOTES:

For high-dimensional systems prefer using AN in sparse format. In this case, the matrix exponential is computed using scipy.sparse.linalg.expm_multiply.

EXAMPLES:

sage: from carlin.transformation import get_Fj_from_model, truncated_matrix
sage: from carlin.library import quadratic_scalar
sage: S = quadratic_scalar()
sage: (f, n, k) = S.funcs(), S.dim(), S.degree()
sage: Fjnk = get_Fj_from_model(f, n, k)

Consider a fourth order truncation:

sage: AN_sparse = truncated_matrix(4, *Fjnk, input_format="Fj_matrices")
sage: AN_sparse.toarray()
array([[ 1.,  1.,  0.,  0.],
       [ 0.,  2.,  2.,  0.],
       [ 0.,  0.,  3.,  3.],
       [ 0.,  0.,  0.,  4.]])

We can solve the linear ODE using the sparse matrix AN_sparse:

sage: from carlin.io import solve_ode_exp
sage: ans = solve_ode_exp(AN_sparse, x0=[0.1], N=4, tini=0, \
            T=1, NPOINTS=20)

It can also be solved using Sage matrices (although this is often less performant, since it works with dense matrices):

sage: AN_dense = matrix(AN_sparse.toarray())
sage: ans = solve_ode_exp(AN_dense, x0=[0.1], N=4, tini=0, \
            T=1, NPOINTS=20)

Polynomial ODE class

A class to represent a system of polynomial ordinary differential equations (ODEs).

AUTHOR:

  • Marcelo Forets (May 2017 at VERIMAG - UGA)
class carlin.polynomial_ode.PolynomialODE(f=None, n=None, k=None)

Bases: sage.structure.sage_object.SageObject

This class represents a finite set of polynomial ODEs.

deg()

Get the degree (max over the degree of all polynomials).

degree()

Get the degree (max over the degree of all polynomials).

dim()

Get the dimesion (number of variables).

funcs()

Get the right-hand side functions list.

plot_solution(x0=None, tini=0, T=1, NPOINTS=100, xcoord=0, ycoord=1, plotjoined=True, **kwargs)

Solve and plot for the given coordinates.

INPUT:

  • x0 – vector; initial condition
  • tini – initial time of simulation
  • T – final time of simulation
  • NPOINTS – number of points sampled
  • xcoord – (default: \(0\)), x-coordinate in plot
  • ycoord – (default: \(1\)), y coordinate in plot

EXAMPLES:

sage: from carlin.library import vanderpol
sage: S = vanderpol(1, 1)
sage: S.plot_solution(x0=[0.5, 1], T=20, NPOINTS=200) # not tested
Graphics object consisting of 1 graphics primitive
solve(x0=None, tini=0, T=1, NPOINTS=100)

Solve the polynomial ODE using GSL.

INPUT:

  • model – PolynomialODE, defining the tuple \((f, n, k)\)
  • N – integer; truncation order
  • x0 – vector; initial condition
  • tini – initial time of simulation
  • T – final time of simulation
  • NPOINTS – number of points sampled
  • xcoord – (default: \(0\)), x-coordinate in plot
  • ycoord – (default: \(1\)), y coordinate in plot

EXAMPLES:

Let us solve the scalar equation \(x' = -2x + 3\):

sage: from carlin.polynomial_ode import PolynomialODE
sage: x = polygen(QQ, 'x')
sage: P = PolynomialODE([-2*x+3], 1, 1)
sage: Solution = P.solve(x0=[0], tini=0, T=4)
sage: type(Solution)
<class 'sage.calculus.ode.ode_solver'>

Let us compute the solution of the vanderpol ODE:

sage: from carlin.library import vanderpol
sage: S = vanderpol(1, 1).solve(x0=[0.5, 1.])
sage: S.solution[1]
(0.01, [0.5100..., 1.0024...])

Library of polynomial ODEs

Library of commonly used or famous polynomial ODE systems.

The following functions are available:

arrowsmith_and_place_fig_3_5e_page_79() Nonlinear two-dimensional system with an hyperbolic fixed point
biomodel_2() Nine-dimensional polynomial ODE form a biological model
chen_seven_dim() Seven-dimensional nonlinear system of quadratic order
cubic_scalar() Scalar ODE with a cubic term
quadratic_scalar() Scalar ODE with a quadratic term
vanderpol() Van der Pol oscillator

AUTHOR:

  • Marcelo Forets (May 2017)
carlin.library.arrowsmith_and_place_fig_3_5e_page_79()

Nonlinear two-dimensional system with an hyperbolic fixed point.

It is defined as:

\[\begin{split}\begin{aligned} x' &= x^2+(x+y)/2 \\ y' &= (-x+3y)/2 \end{aligned}\end{split}\]

Taken from p. 79 of the book by Arrowsmith and Place, Dynamical Systems: Differential Equations, maps and chaotic behaviour.

carlin.library.biomodel_2()

This is a nine-dimensional polynomial ODE used as benchmark model in the Flow star tool.

The model is adapted from E. Klipp, R. Herwig, A. Kowald, C. Wierling, H. Lehrach. Systems Biology in Practice: Concepts, Implementation and Application. Wiley-Blackwell, 2005.

carlin.library.chen_seven_dim(u=0)

This is a seven-dimensional nonlinear system of quadratic order.

It appears as 'example_nonlinear_reach_04_sevenDim_nonConvexRepr.m' in the tool CORA 2016, in the examples for continuous nonlinear systems.

There is an independent term, \(u\), in the fourth equation with value \(2.0\) that has been neglected here for convenience (hence we take \(u=0\) by default).

carlin.library.cubic_scalar(a=1, b=1)

A scalar ODE with a cubic term.

It is defined as:

\[x'(t) = ax(t) + bx(t)^3\]

where \(a\) and \(b\) are paremeters of the ODE.

EXAMPLES:

sage: from carlin.library import cubic_scalar
sage: C = cubic_scalar(-1, 1)
sage: C.funcs()
[x0^3 - x0]

Compute the Carleman embedding truncated at order \(N=4\):

sage: from carlin.transformation import get_Fj_from_model, truncated_matrix 
sage: Fj = get_Fj_from_model(C.funcs(), C.dim(), C.degree())
sage: matrix(truncated_matrix(4, *Fj, input_format="Fj_matrices").toarray())
[-1.0  0.0  1.0  0.0]
[ 0.0 -2.0  0.0  2.0]
[ 0.0  0.0 -3.0  0.0]
[ 0.0  0.0  0.0 -4.0]
carlin.library.quadratic_scalar(a=1, b=1)

A scalar ODE with a quadratic term.

It is defined as:

\[x'(t) = ax(t) + bx(t)^2\]

where \(a\) and \(b\) are paremeters of the ODE.

EXAMPLES:

sage: from carlin.library import quadratic_scalar
sage: Q = quadratic_scalar(); Q
A Polynomial ODE in n = 1 variables
sage: Q.funcs()
[x0^2 + x0]

Compute the Carleman embedding truncated at order \(N=4\):

sage: from carlin.transformation import get_Fj_from_model, truncated_matrix 
sage: Fj = get_Fj_from_model(Q.funcs(), Q.dim(), Q.degree())
sage: matrix(truncated_matrix(4, *Fj, input_format="Fj_matrices").toarray())
[1.0 1.0 0.0 0.0]
[0.0 2.0 2.0 0.0]
[0.0 0.0 3.0 3.0]
[0.0 0.0 0.0 4.0]
carlin.library.vanderpol(mu=1, omega=1)

The Van der Pol oscillator is a non-conservative system with non-linear damping.

It is defined as:

\[\begin{split}\begin{aligned} x' &= y \\ y' &= -\omega^2 x - (x^2 - 1) \mu y \end{aligned}\end{split}\]

where \(\omega\) is the natural frequency and \(\mu\) is the damping parameter. For additional information see the Wikipedia article Van_der_Pol_oscillator.

EXAMPLES:

sage: from carlin.library import vanderpol
sage: vanderpol(SR.var('mu'), SR.var('omega')).funcs()
[x1, -omega^2*x0 - (x0^2 - 1)*mu*x1]

Polyhedron Tools

See polyhedron_tools