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 thetoarray
method is availablen
– dimension on state-spacek
– order of the systemord
– order of the \(p\)-th norm, for \(1 \leq p < \infty\), andp='inf'
for \(p=\infty\)
OUTPUT:
Dictionary
c
containingnorm_Fi_inf
,log_norm_F1_inf
andbeta0_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 formatN
– integer; truncation orderx0
– list; initial point
OUTPUT:
Ts
– convergence time computed from the reduced quadratic systemerror
– 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 formattarget_filename
– string with the name of the output file in MAT formatN
– truncation orderx0
– 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 formatN
– truncation orderx0
– initial point, a listTfrac
– (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’sn
– integer, system’s dimensionk
– 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 truncationF
– list, sequence of matrices \(F_j\)n
– the dimension of the state-spacek
– the order of the polynomial vector field. It is equal tolen(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 truncationinput_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 pathF
– list of sparse matrices \(F_1, \ldots, F_k\)n
– dimension of state-spacek
– 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 systemk
– 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 fk
– 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 orderx0
– vector; initial conditiontini
– initial time of simulationT
– final time of simulationNPOINTS
– number of points sampledxcoord
– (default: \(0\)), x-coordinate in plotycoord
– (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 formatx0
– vector, initial pointN
– integer, order of the truncationtini
– (optional, default: 0) initial timeT
– (optional, default: 1) final timeNPOINTS
– (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 conditiontini
– initial time of simulationT
– final time of simulationNPOINTS
– number of points sampledxcoord
– (default: \(0\)), x-coordinate in plotycoord
– (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 orderx0
– vector; initial conditiontini
– initial time of simulationT
– final time of simulationNPOINTS
– number of points sampledxcoord
– (default: \(0\)), x-coordinate in plotycoord
– (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