In mathematics and physics, the Magnus expansion, named after Wilhelm Magnus (1907–1990), provides an exponential representation of the product integral solution of a first-order homogeneous linear differential equation for a linear operator. In particular, it furnishes the fundamental matrix of a system of linear ordinary differential equations of order n with varying coefficients. The exponent is aggregated as an infinite series, whose terms involve multiple integrals and nested commutators.
The deterministic case
Magnus approach and its interpretation Given the n × n coefficient matrix A(t), one wishes to solve the initial-value problem associated with the linear ordinary differential equation
Y ′ ( t ) = A ( t ) Y ( t ) , Y ( t 0 ) = Y 0 {\displaystyle Y'(t)=A(t)Y(t),\quad Y(t_{0})=Y_{0}}
for the unknown n-dimensional vector function Y(t). When n = 1, the solution is given as a product integral
Y ( t ) = exp ( ∫ t 0 t A ( s ) d s ) Y 0 . {\displaystyle Y(t)=\exp \left(\int _{t_{0}}^{t}A(s)\,ds\right)Y_{0}.}
This is still valid for n > 1 if the matrix A(t) satisfies A(t1) A(t2) = A(t2) A(t1) for any pair of values of t, t1 and t2. In particular, this is the case if the matrix A is independent of t. In the general case, however, the expression above is no longer the solution of the problem. The approach introduced by Magnus to solve the matrix initial-value problem is to express the solution by means of the exponential of a certain n × n matrix function Ω(t, t0):
Y ( t ) = exp ( Ω ( t , t 0 ) ) Y 0 , {\displaystyle Y(t)=\exp {\big (}\Omega (t,t_{0}){\big )}\,Y_{0},}
which is subsequently constructed as a series expansion:
Ω ( t ) = ∑ k = 1 ∞ Ω k ( t ) , {\displaystyle \Omega (t)=\sum _{k=1}^{\infty }\Omega _{k}(t),}
where, for simplicity, it is customary to write Ω(t) for Ω(t, t0) and to take t0 = 0. Magnus appreciated that, since d/dt (eΩ) e−Ω = A(t), using a Poincaré−Hausdorff matrix identity, he could relate the time derivative of Ω to the generating function of Bernoulli numbers and the adjoint endomorphism of Ω,
Ω ′ = ad Ω exp ( ad Ω ) − 1 A , {\displaystyle \Omega '={\frac {\operatorname {ad} _{\Omega }}{\exp(\operatorname {ad} _{\Omega })-1}}A,}
to solve for Ω recursively in terms of A "in a continuous analog of the BCH expansion", as outlined in a subsequent section. The equation above constitutes the Magnus expansion, or Magnus series, for the solution of matrix linear initial-value problem. The first four terms of this series read
Ω 1 ( t ) = ∫ 0 t A ( t 1 ) d t 1 , Ω 2 ( t ) = 1 2 ∫ 0 t d t 1 ∫ 0 t 1 d t 2 [ A ( t 1 ) , A ( t 2 ) ] , Ω 3 ( t ) = 1 6 ∫ 0 t d t 1 ∫ 0 t 1 d t 2 ∫ 0 t 2 d t 3 ( [ A ( t 1 ) , [ A ( t 2 ) , A ( t 3 ) ] ] + [ A ( t 3 ) , [ A ( t 2 ) , A ( t 1 ) ] ] ) , Ω 4 ( t ) = 1 12 ∫ 0 t d t 1 ∫ 0 t 1 d t 2 ∫ 0 t 2 d t 3 ∫ 0 t 3 d t 4 ( [ [ [ A 1 , A 2 ] , A 3 ] , A 4 ] + [ A 1 , [ [ A 2 , A 3 ] , A 4 ] ] + [ A 1 , [ A 2 , [ A 3 , A 4 ] ] ] + [ A 2 , [ A 3 , [ A 4 , A 1 ] ] ] ) , {\displaystyle {\begin{aligned}\Omega _{1}(t)&=\int _{0}^{t}A(t_{1})\,dt_{1},\\\Omega _{2}(t)&={\frac {1}{2}}\int _{0}^{t}dt_{1}\int _{0}^{t_{1}}dt_{2}\,[A(t_{1}),A(t_{2})],\\\Omega _{3}(t)&={\frac {1}{6}}\int _{0}^{t}dt_{1}\int _{0}^{t_{1}}dt_{2}\int _{0}^{t_{2}}dt_{3}\,{\Bigl (}{\big [}A(t_{1}),[A(t_{2}),A(t_{3})]{\big ]}+{\big [}A(t_{3}),[A(t_{2}),A(t_{1})]{\big ]}{\Bigr )},\\\Omega _{4}(t)&={\frac {1}{12}}\int _{0}^{t}dt_{1}\int _{0}^{t_{1}}dt_{2}\int _{0}^{t_{2}}dt_{3}\int _{0}^{t_{3}}dt_{4}\,\left({\Big [}{\big [}[A_{1},A_{2}],A_{3}{\big ]},A_{4}{\Big ]}\right.\\&\qquad +{\Big [}A_{1},{\big [}[A_{2},A_{3}],A_{4}{\big ]}{\Big ]}+{\Big [}A_{1},{\big [}A_{2},[A_{3},A_{4}]{\big ]}{\Big ]}+\left.{\Big [}A_{2},{\big [}A_{3},[A_{4},A_{1}]{\big ]}{\Big ]}\right),\end{aligned}}}
where [A, B] ≡ A B − B A is the matrix commutator of A and B. These equations may be interpreted as follows: Ω1(t) coincides exactly with the exponent in the scalar (n = 1) case, but this equation cannot give the whole solution. If one insists in having an exponential representation (Lie group), the exponent needs to be corrected. The rest of the Magnus series provides that correction systematically: Ω or parts of it are in the Lie algebra of the Lie group on the solution. In applications, one can rarely sum exactly the Magnus series, and one has to truncate it to get approximate solutions. The main advantage of the Magnus proposal is that the truncated series very often shares important qualitative properties with the exact solution, similarly to other conventional perturbation theories. For instance, in classical mechanics the symplectic character of the time evolution is preserved at every order of approximation. Similarly, the unitary character of the time evolution operator in quantum mechanics is also preserved (in contrast, e.g., to the Dyson series solving the same problem).
Convergence of the expansion From a mathematical point of view, the convergence problem is the following: given a certain matrix A(t), when can the exponent Ω(t) be obtained as the sum of the Magnus series? A sufficient condition for this series to converge for t ∈ [0,T) is
∫ 0 T ‖ A ( s ) ‖ 2 d s < π , {\displaystyle \int _{0}^{T}\|A(s)\|_{2}\,ds<\pi ,}
where ‖ ⋅ ‖ 2 {\displaystyle \|\cdot \|_{2}} denotes a matrix norm. This result is generic in the sense that one may construct specific matrices A(t) for which the series diverges for any t > T.
Magnus generator A recursive procedure to generate all the terms in the Magnus expansion utilizes the matrices Sn(k) defined recursively through
S n ( j ) = ∑ m = 1 n − j [ Ω m , S n − m ( j − 1 ) ] , 2 ≤ j ≤ n − 1 , {\displaystyle S_{n}^{(j)}=\sum _{m=1}^{n-j}\left[\Omega _{m},S_{n-m}^{(j-1)}\right],\quad 2\leq j\leq n-1,}
S n ( 1 ) = [ Ω n − 1 , A ] , S n ( n − 1 ) = ad Ω 1 n − 1 ( A ) , {\displaystyle S_{n}^{(1)}=\left[\Omega _{n-1},A\right],\quad S_{n}^{(n-1)}=\operatorname {ad} _{\Omega _{1}}^{n-1}(A),}
which then furnish
Ω 1 = ∫ 0 t A ( τ ) d τ , {\displaystyle \Omega _{1}=\int _{0}^{t}A(\tau )\,d\tau ,}
Ω n = ∑ j = 1 n − 1 B j j ! ∫ 0 t S n ( j ) ( τ ) d τ , n ≥ 2. {\displaystyle \Omega _{n}=\sum _{j=1}^{n-1}{\frac {B_{j}}{j!}}\int _{0}^{t}S_{n}^{(j)}(\tau )\,d\tau ,\quad n\geq 2.}
Here adkΩ is a shorthand for an iterated commutator (see adjoint endomorphism):
ad Ω 0 A = A , ad Ω k + 1 A = [ Ω , ad Ω k A ] , {\displaystyle \operatorname {ad} _{\Omega }^{0}A=A,\quad \operatorname {ad} _{\Omega }^{k+1}A=[\Omega ,\operatorname {ad} _{\Omega }^{k}A],}
while Bj are the Bernoulli numbers with B1 = −1/2. Finally, when this recursion is worked out explicitly, it is possible to express Ωn(t) as a linear combination of n-fold integrals of n − 1 nested commutators involving n matrices A:
Ω n ( t ) = ∑ j = 1 n − 1 B j j ! ∑ k 1 + ⋯ + k j = n − 1 k 1 ≥ 1 , … , k j ≥ 1 ∫ 0 t ad Ω k 1 ( τ ) ad Ω k 2 ( τ ) ⋯ ad Ω k j ( τ ) A ( τ ) d τ , n ≥ 2 , {\displaystyle \Omega _{n}(t)=\sum _{j=1}^{n-1}{\frac {B_{j}}{j!}}\sum _{k_{1}+\cdots +k_{j}=n-1 \atop k_{1}\geq 1,\ldots ,k_{j}\geq 1}\int _{0}^{t}\operatorname {ad} _{\Omega _{k_{1}}(\tau )}\operatorname {ad} _{\Omega _{k_{2}}(\tau )}\cdots \operatorname {ad} _{\Omega _{k_{j}}(\tau )}A(\tau )\,d\tau ,\quad n\geq 2,}
which becomes increasingly intricate with n.
Canonical Representation A non-recursive expression for all the terms in the Magnus expansion can also be formulated by encoding the combinatorics of the recursion in a single integral. This was found by Bialynicki-Birula, Mielnik, and Plebañski in 1969. The iterative solution Y(t,t_0) to the initial value problem can be written as: Y ( t , t 0 ) = I + ∑ n = 1 ∞ ∫ t 0 t d t n … ∫ t 0 t A ( t n ) θ n , n − 1 d t 1 A ( t n − 1 ) … θ 2 , 1 A ( t 1 ) θ k , l = { 1 , t k ≥ t l 0 , t k < t l {\displaystyle Y(t,t_{0})=I+\sum _{n=1}^{\infty }\int _{t_{0}}^{t}dt_{n}\dots \int _{t_{0}}^{t}A(t_{n})\theta _{n,n-1}dt_{1}A(t_{n-1})\dots \theta _{2,1}A(t_{1})\qquad \theta _{k,l}={\begin{cases}1,\quad t_{k}\geq t_{l}\\0,\quad t_{k}<t_{l}\end{cases}}} Then, this representation can be used to derive the canonical representation of Ω n {\displaystyle \Omega _{n}} , using a function Θ n ( t 1 , … , t n ) {\displaystyle \Theta _{n}(t_{1},\dots ,t_{n})} which counts the time ordering, and thus encodes the combinatorial structure of the expansion:
Θ n ( t 1 , t 2 , … , t n ) = θ n , n − 1 + θ n − 1 , n − 2 + … θ 2 , 1 {\displaystyle \Theta _{n}(t_{1},t_{2},\dots ,t_{n})=\theta _{n,n-1}+\theta _{n-1,n-2}+\dots \theta _{2,1}}
Ω n = 1 n ! ∫ t 0 t d t n … ∫ t 0 t d t 1 ( − 1 ) n − Θ − 1 Θ ! ( n − Θ − 1 ) ! A ( t n ) ⋯ A ( t 1 ) {\displaystyle \Omega _{n}={\frac {1}{n!}}\int _{t_{0}}^{t}dt_{n}\dots \int _{t_{0}}^{t}dt_{1}(-1)^{n-\Theta -1}\Theta !(n-\Theta -1)!A(t_{n})\cdots A(t_{1})}
This format is advantageous in many applications due to its independent limits of integration and closed form for an arbitrary degree.
The stochastic case
Extension to stochastic ordinary differential equations For the extension to the stochastic case let ( W t ) t ∈ [ 0 , T ] {\textstyle \left(W_{t}\right)_{t\in [0,T]}} be a R q {\textstyle \mathbb {R} ^{q}} -dimensional Brownian motion, q ∈ N > 0 {\textstyle q\in \mathbb {N} _{>0}} , on the probability space ( Ω , F , P ) {\textstyle \left(\Omega ,{\mathcal {F}},\mathbb {P} \right)}
with finite time horizon T > 0 {\textstyle T>0} and natural filtration. Now, consider the linear matrix-valued stochastic Itô differential equation (with Einstein's summation convention over the index j)
d X t = B t X t d t + A t ( j ) X t d W t j , X 0 = I d , d ∈ N > 0 , {\displaystyle dX_{t}=B_{t}X_{t}dt+A_{t}^{(j)}X_{t}dW_{t}^{j},\quad X_{0}=I_{d},\qquad d\in \mathbb {N} _{>0},}
where B ⋅ , A ⋅ ( 1 ) , … , A ⋅ ( j ) {\textstyle B_{\cdot },A_{\cdot }^{(1)},\dots ,A_{\cdot }^{(j)}} are progressively measurable d × d {\textstyle d\times d} -valued bounded stochastic processes and I d {\textstyle I_{d}} is the identity matrix. Following the same approach as in the deterministic case with alterations due to the stochastic setting the corresponding matrix logarithm will turn out as an Itô-process, whose first two expansion orders are given by Y t ( 1 ) = Y t ( 1 , 0 ) + Y t ( 0 , 1 ) {\textstyle Y_{t}^{(1)}=Y_{t}^{(1,0)}+Y_{t}^{(0,1)}} and Y t ( 2 ) = Y t ( 2 , 0 ) + Y t ( 1 , 1 ) + Y t ( 0 , 2 ) {\textstyle Y_{t}^{(2)}=Y_{t}^{(2,0)}+Y_{t}^{(1,1)}+Y_{t}^{(0,2)}} , where with Einstein's summation convention over i and j
Y t ( 0 , 0 ) = 0 , Y t ( 1 , 0 ) = ∫ 0 t A s ( j ) d W s j , Y t ( 0 , 1 ) = ∫ 0 t B s d s , Y t ( 2 , 0 )