In mathematics, Carleman linearization (also called Carleman embedding) is a method that formally represents a finite-dimensional nonlinear system of ordinary differential equations as an infinite-dimensional linear system. Truncating the resulting hierarchy at a finite order yields approximate solutions of the original nonlinear system, known as Carleman approximants. The method was introduced by the Swedish mathematician Torsten Carleman in 1932. For polynomial systems, the construction successively introduces monomials of the state variables as new dependent variables. More generally, the additional variables may be chosen according to the nonlinear structure of the differential system. Carleman linearization may be regarded as a systematic extension of the usual local linearization of nonlinear differential equations, in which higher-order nonlinear terms are retained through an enlargement of the phase space. Carleman linearization is related to composition operators, particularly the Koopman operator, which represents nonlinear dynamics through linear evolution on a space of observable functions. It has been applied in nonlinear control theory, state estimation, chemical kinetics, process modelling, mathematical biology, and the construction of finite-dimensional approximations to Koopman operators. The method has also been used in the development of quantum algorithms for nonlinear differential equations.
Method Consider an autonomous system of ordinary differential equations,
d x d t = f ( x ) , x ( 0 ) = c , x ∈ R p . {\displaystyle {\frac {d\mathbf {x} }{dt}}=\mathbf {f} (\mathbf {x} ),\qquad \mathbf {x} (0)=\mathbf {c} ,\qquad \mathbf {x} \in \mathbb {R} ^{p}.}
Carleman linearization consists of replacing nonlinear terms appearing on the right-hand side of the system with new dependent variables. Differentiating these variables generally produces further nonlinear terms, which are treated successively in the same way. The procedure therefore generates an infinite sequence of linear differential equations. A finite-dimensional approximation is obtained when the sequence is terminated according to a chosen truncation criterion. The solution of the truncated linear system is called a Carleman approximant. For polynomial systems, the additional variables are monomials of the original state variables. The order N {\displaystyle N} of the approximation is then defined as the degree of the highest-order monomial retained. Other sets of dependent variables may be chosen when they are better adapted to the structure of the differential system. After truncation at order N {\displaystyle N} , the original p {\displaystyle p} -dimensional nonlinear system is replaced by a finite-dimensional linear system of dimension m ≥ p {\displaystyle m\geq p} . The value of m {\displaystyle m} depends on the approximation order and on the nonlinear structure of the original equations. The approximation to the original solution is given by the components of the enlarged linear system corresponding to the original state variables. Two forms of the construction can be distinguished according to the point in phase space around which the Carleman approximant is built: linearization around an equilibrium point, referred to as Type A, and linearization around the initial condition, referred to as Type B.
Type A: linearization around an equilibrium point Let q {\displaystyle \mathbf {q} } be an equilibrium point of the nonlinear system,
f ( q ) = 0 . {\displaystyle \mathbf {f} (\mathbf {q} )=\mathbf {0} .}
Introducing the shifted variable
y = x − q {\displaystyle \mathbf {y} =\mathbf {x} -\mathbf {q} }
places the equilibrium point at the origin. Since the constant term of the transformed differential system vanishes, truncation of the Carleman sequence produces a homogeneous linear system with constant coefficients,
d X N d t = M N X N , {\displaystyle {\frac {d\mathbf {X} _{N}}{dt}}=M_{N}\mathbf {X} _{N},}
where X N {\displaystyle \mathbf {X} _{N}} contains the retained variables and M N {\displaystyle M_{N}} is the corresponding m × m {\displaystyle m\times m} coefficient matrix. Its solution is
X N ( t ) = e M N t X N ( 0 ) . {\displaystyle \mathbf {X} _{N}(t)=e^{M_{N}t}\mathbf {X} _{N}(0).}
The approximate solution of the original nonlinear system is provided by the p {\displaystyle p} components of X N ( t ) {\displaystyle \mathbf {X} _{N}(t)} associated with the original variables, after reversing the coordinate shift. The first-order Type A approximant coincides with the usual local linearization obtained from the Jacobian matrix of f {\displaystyle \mathbf {f} } at the equilibrium point. Higher-order approximants include additional nonlinear terms and may therefore be regarded as an extension of the local or qualitative analysis of nonlinear differential systems. For a fixed approximation order, their accuracy generally depends on the distance between the trajectory of interest and the equilibrium point around which the approximant was constructed.
Type B: linearization around the initial condition Alternatively, the Carleman approximant may be constructed around the initial condition itself. Introducing
y = x − c , y ( 0 ) = 0 , {\displaystyle \mathbf {y} =\mathbf {x} -\mathbf {c} ,\qquad \mathbf {y} (0)=\mathbf {0} ,}
gives
d y d t = f ( c + y ) . {\displaystyle {\frac {d\mathbf {y} }{dt}}=\mathbf {f} (\mathbf {c} +\mathbf {y} ).}
Unless c {\displaystyle \mathbf {c} } is an equilibrium point, the transformed system contains the nonzero constant term f ( c ) {\displaystyle \mathbf {f} (\mathbf {c} )} . Consequently, truncation of the Carleman sequence produces a non-homogeneous linear system,
d Y N d t = M N Y N + z N , Y N ( 0 ) = 0 , {\displaystyle {\frac {d\mathbf {Y} _{N}}{dt}}=M_{N}\mathbf {Y} _{N}+\mathbf {z} _{N},\qquad \mathbf {Y} _{N}(0)=\mathbf {0} ,}
where z N {\displaystyle \mathbf {z} _{N}} is a constant vector determined by the initial condition and by the nonlinear differential system. The solution can be written as
Y N ( t ) = ∫ 0 t e M N ( t − s ) z N d s . {\displaystyle \mathbf {Y} _{N}(t)=\int _{0}^{t}e^{M_{N}(t-s)}\mathbf {z} _{N}\,ds.}
If M N {\displaystyle M_{N}} is invertible, this becomes
Y N ( t ) = M N − 1 ( e M N t − I ) z N . {\displaystyle \mathbf {Y} _{N}(t)=M_{N}^{-1}\left(e^{M_{N}t}-I\right)\mathbf {z} _{N}.}
The approximation to x ( t ) {\displaystyle \mathbf {x} (t)} is obtained by adding c {\displaystyle \mathbf {c} } to the components of Y N ( t ) {\displaystyle \mathbf {Y} _{N}(t)} associated with the shifted original variables. Type B approximants require the solution of a non-homogeneous linear system and generally involve a greater algebraic burden than Type A approximants. However, centring the approximation at the initial condition can substantially improve its local accuracy.
Local accuracy Comparison of the derivatives at t = 0 {\displaystyle t=0} of the exact solution and the Carleman approximants reveals different local-accuracy patterns. In the one- and two-dimensional systems analysed, an approximant of order N {\displaystyle N} reproduces the derivatives of the exact solution from order 0 {\displaystyle 0} through order N − 1 {\displaystyle N-1} , independently of the point around which it is constructed. In one dimension, let q {\displaystyle q} be the point around which the approximant x ~ ( t ) {\displaystyle {\widetilde {x}}(t)} is constructed, and let c = x ( 0 ) {\displaystyle c=x(0)} be the initial condition. The differences between the derivatives of the exact and approximate solutions have the following pattern:
When the approximant is centred at the initial condition, q = c {\displaystyle q=c} , the differences in the second row also vanish. Thus, for the one-dimensional systems considered, a Type B approximant reproduces the derivatives of the exact solution through order 2 N {\displaystyle 2N} . In two dimensions, let the approximants x ~ 1 ( t ) {\displaystyle {\widetilde {x}}_{1}(t)} and x ~ 2 ( t ) {\displaystyle {\widetilde {x}}_{2}(t)} be constructed around ( q 1 , q 2 ) {\displaystyle (q_{1},q_{2})} , and let ( c 1 , c 2 ) {\displaystyle (c_{1},c_{2})} be the initial condition. For i = 1 , 2 {\displaystyle i=1,2} , the corresponding pattern is
When ( q 1 , q 2 ) = ( c 1 , c 2 ) {\displaystyle (q_{1},q_{2})=(c_{1},c_{2})} , all the displayed differences vanish. The additional agreement of the derivatives from order N {\displaystyle N} through order 2 N {\displaystyle 2N} accounts for the greater local accuracy observed for Type B approximants of a given order.
Examples The following examples illustrate the construction of Type A and Type B Carleman approximants in one and two dimensions. Only one equilibrium point is considered in each Type A construction.
One-dimensional example: Riccati equation Consider the Riccati equation
u ˙ = − u + u 2 , u ( 0 ) = c . {\displaystyle {\dot {u}}=-u+u^{2},\qquad u(0)=c.}
Its exact solution is
u ( t ) = c c + ( 1 − c ) e t . {\displaystyle u(t)={\frac {c}{c+(1-c)e^{t}}}.}
The equation has two equilibrium points, u = 0 {\displaystyle u=0} and u = 1 {\displaystyle u=1} . The Type A approximant below is constructed only around u = 0 {\displaystyle u=0} .
Type A approximant Introduce the Carleman variables
u k ( t ) = u k ( t ) , k = 1 , 2 , … . {\displaystyle u_{k}(t)=u^{k}(t),\qquad k=1,2,\ldots .}
Differentiation gives the infinite linear sequence
u ˙ k = − k u k + k u k + 1 , k = 1 , 2 , … . {\displaystyle {\dot {u}}_{k}=-ku_{k}+ku_{k+1},\qquad k=1,2,\ldots .}
A second-order approximation is obtained by retaining u 1 = u {\displaystyle u_{1}=u} and u 2 = u 2 {\displaystyle u_{2}=u^{2}} and discarding the term u 3 {\displaystyle u_{3}} . The resulting homogeneous system is
d d t ( u 1 u 2 ) = ( − 1 1 0 − 2 ) ( u 1 u 2 ) , ( u 1 ( 0 ) u 2 ( 0 ) ) = ( c c 2 ) . {\displaystyle {\frac {d}{dt}}{\begin{pmatrix}u_{1}\\u_{2}\end{pmatrix}}={\begin{pmatrix}-1&1\\0&-2\end{pmatrix}}{\begin{pmatrix}u_{1}\\u_{2}\end{pmatrix}},\qquad {\begin{pmatrix}u_{1}(0)\\u_{2}(0)\end{pmatrix}}={\begin{pmatrix}c\\c^{2}\end{pmatrix}}.}
Solving the linear system gives the second-order Type A approximant
u ~ N = 2 A ( t ) = u 1 ( t ) = c [ 1 + c ( 1 − e − t ) ] e − t . {\displaystyle {\widetilde {u}}_{N=2}^{\,A}(t)=u_{1}(t)=c\left[1+c\left(1-e^{-t}\right)\right]e^{-t}.}
Type B approximant To construct the approximant around the initial condition, introduce
u = w + c , w ( 0 ) = 0. {\displaystyle u=w+c,\qquad w(0)=0.}
The differential equation becomes
w ˙ = ( 2 c − 1 ) w + w 2 + c ( c − 1 ) . {\displaystyle {\dot {w}}=(2c-1)w+w^{2}+c(c-1).}
Define
α = 2 c − 1 , β = c ( c − 1 ) . {\displaystyle \alpha =2c-1,\qquad \beta =c(c-1).}
Retaining w {\displaystyle w} and w 2 {\displaystyle w^{2}} gives the second-order non-homogeneous linear system
d d t ( w w 2 ) = ( α 1 2 β 2 α ) ( w w 2 ) + ( β 0 ) , ( w ( 0 ) w 2 ( 0 ) ) = ( 0 0 ) . {\displaystyle {\frac {d}{dt}}{\begin{pmatrix}w\\w^{2}\end{pmatrix}}={\begin{pmatrix}\alpha &1\\2\beta &2\alpha \end{pmatrix}}{\begin{pmatrix}w\\w^{2}\end{pmatrix}}+{\begin{pmatrix}\beta \\0\end{pmatrix}},\qquad {\begin{pmatrix}w(0)\\w^{2}(0)\end{pmatrix}}={\begin{pmatrix}0\\0\end{pmatrix}}.}
Writing
M = ( α 1 2 β 2 α ) , b = ( β 0 ) , {\displaystyle M={\begin{pmatrix}\alpha &1\\2\beta &2\alpha \end{pmatrix}},\qquad \mathbf {b} ={\begin{pmatrix}\beta \\0\end{pmatrix}},}
the solution is
( w ( t ) w 2 ( t ) ) = M − 1 ( e M t − I ) b , {\displaystyle {\begin{pmatrix}w(t)\\w^{2}(t)\end{pmatrix}}=M^{-1}\left(e^{Mt}-I\right)\mathbf {b} ,}
when M {\displaystyle M} is invertible. The second-order Type B approximant is therefore
u ~ N = 2 B ( t ) = c + w ( t ) = c { 1 + c − 1 3 c ( c − 1 ) + 1 [ 1 − 2 c + 1 2 ξ ( [ ( 2 c − 1 ) + 6 c ( c − 1 ) + 1 ] e λ − t + [ ( 2 c − 1 ) − 6 c ( c − 1 ) − 1 ] e λ + t ) ] } {\displaystyle {\begin{array}{l}{\widetilde {u}}_{N=2}^{\,B}(t)=c+w(t)=\\[4pt]\displaystyle c{\Biggl \{}1+{\frac {c-1}{3c(c-1)+1}}{\Biggl [}1-2c+{\frac {1}{2\xi }}{\Biggl (}\left[(2c-1)+6c(c-1)+1\right]e^{\lambda _{-}t}\\[4pt]+\left[(2c-1)-6c(c-1)-1\right]e^{\lambda _{+}t}{\Biggr )}{\Biggr ]}{\Biggr \}}\end{array}}}
where ξ ≡ 1 + 12 c ( c − 1 ) {\displaystyle \xi \equiv {\sqrt {1+12c(c-1)}}} , and λ ± = [ 3 ( 2 c − 1 ) ± ξ ] / 2 {\displaystyle \lambda _{\pm }=[3(2c-1)\pm \xi ]/2} are the eigenvalues of M {\displaystyle M} . Although this expression is algebraically more involved than the Type A approximant, centring the construction at u ( 0 ) = c {\displaystyle u(0)=c} generally improves its local accuracy. The plot compares the exact solutions with second- and third-order Type A approximants constructed around the two equilibrium points and with a second-order Type B approximant constructed around each initial condition. [[File:Carleman approximants for the Riccati equation.svg
|thumb |right |upright=1.45 |alt=Riccati equation. Comparison of the exact solutions with second- and third-order Carleman approximants constructed around the equilibrium points, and second-order approximants constructed around the four initial conditions |Riccati equation. Exact solutions (black dotted curves) and Carleman approximants (solid curves) as functions of t {\displaystyle t} for four initial conditions, c = 1.1 , 0.9 , 0.5 , − 0.8 {\displaystyle c=1.1,0.9,0.5,-0.8} . The subscript in the legend indicates the point around which the approximant is centred, and the superscript gives its order. Dashed lines are the two equilibrium solutions.
]]
Two-dimensional example: center Consider the two-dimensional nonlinear system
x ˙ = − y − x y , y ˙ = x + x y , {\displaystyle {\dot {x}}=-y-xy,\qquad {\dot {y}}=x+xy,}
with initial condition
x ( 0 ) = a , y ( 0 ) = b . {\displaystyle x(0)=a,\qquad y(0)=b.}
The system has an equilibrium at ( 0 , 0 ) {\displaystyle (0,0)} , which is a center and is stable in the sense of Lyapunov, but not asymptotically stable. It also has a saddle point at ( − 1 , − 1 ) {\displaystyle (-1,-1)} . For x , y > − 1 {\displaystyle x,y>-1} , its trajectories satisfy the first integral
x − ln ( 1 + x ) + y − ln ( 1 + y ) = C , {\displaystyle x-\ln(1+x)+y-\ln(1+y)=C,}
where C {\displaystyle C} is constant along each trajectory.
Type A approximant For the Type A construction, consider only the equilibrium point ( 0 , 0 ) {\displaystyle (0,0)} . At second order, introduce the enlarged set of variables
X = ( y y 2 x x y x 2 ) . {\displaystyle \mathbf {X} ={\begin{pmatrix}y\\y^{2}\\x\\xy\\x^{2}\end{pmatrix}}.}
After discarding terms of degree greater than two, the Carleman system is
d d t ( y y 2 x x y x 2 ) = ( 0 0 1 1 0 0 0 0 2 0 − 1