Continuous changes of variables and the Magnus expansion

In this paper, we are concerned with a formulation of Magnus and Floquet-Magnus expansions for general nonlinear differential equations. To this aim, we introduce suitable continuous variable transformations generated by operators. As an application of the simple formulas so-obtained, we explicitly compute the first terms of the Floquet-Magnus expansion for the Van der Pol oscillator and the nonlinear Schr\"odinger equation on the torus.


Introduction
The Magnus expansion constitutes nowadays a standard tool for obtaining both analytic and numerical approximations to the solutions of non-autonomous linear differential equations. In its simplest formulation, the Magnus expansion [28] aims to construct the solution of the linear differential equation where A(t) is a n × n matrix, as where Ω is an infinite series with whose terms are increasingly complex expressions involving iterated integrals of nested commutators of the matrix A evaluated at different times.
Since the 1960s the Magnus expansion (often with different names) has been used in many different fields, ranging from nuclear, atomic and molecular physics to nuclear magnetic resonance and quantum electrodynamics, mainly in connection with perturbation theory. More recently, it has also been the starting point to construct numerical integration methods in the realm of geometric numerical integration (see [7] for a review), when preserving the main qualitative features of the exact solution, such as its invariant quantities or the geometric structure is at issue [5,23]. The convergence of the expansion is also an important feature and several general results are available [6,10,31,27].
Given the favourable properties exhibited by the Magnus expansion in the treatment of the linear problem (1), it comes as no surprise that several generalizations have been proposed along the years. We can mention, in particular, equation (1) when the (in general complex) matrix-valued function A(t) is periodic with period T . In that case, it is possible to combine the Magnus expansion with the Floquet theorem [16] and construct the solution as where Λ(t + T ) = Λ(t) and both Λ(t) and F are series expansions with Λ k (0) = 0 for all k. This is the so-called Floquet-Magnus expansion [11], and has been widely used in problems of solid state physics and nuclear magnetic resonance [26,29]. Notice that, due to the periodicity of Λ k , the constant term F n can be independently obtained as F k = Ω k (T )/T for all k.
In the general case of a nonlinear ordinary differential equation in R n , the usual procedure to construct the Magnus expansion requires first to transform (6) into a certain linear equation involving operators [1]. This is done by introducing the Lie derivative associated with g and the family of linear transformations Φ t such that where ϕ t denotes the exact flow defined by (6) and f is any (infinitely) differentiable map f : R n −→ R. The operator Φ t obeys a linear differential equation which is then formally solved with the corresponding Magnus expansion [7]. Once the series is truncated, it corresponds to the Lie derivative of some function W (x, t).
Finally, the solution at some given time t = T can be approximated by determining the 1-flow of the autonomous differential equation since, by construction, y(1) ϕ T (x 0 ). Clearly, the whole procedure is different and more involved than in the linear case. It is the purpose of this work to provide a unified framework to derive the Magnus expansion in a simpler way without requiring the apparatus of chronological calculus. This will be possible by applying the continuous transformation theory developed by Dewar in perturbation theory in classical mechanics [17]. In that context, the Magnus series is just the generator of the continuous transformation sending the original system (6) to the trivial one X = 0. Moreover, the same idea can be applied to the Floquet-Magnus expansion, thus establishing a natural connection with the stroboscopic averaging formalism. In the process, the relation with pre-Lie algebras and other combinatorial objects will appear in a natural way. The plan of the paper is as follows. We review several procedures to derive the Magnus expansion for the linear equation (1) in section 2 and introduce a binary operator that will play an important role in the sequel. In section 3 we consider continuous changes of variables and their generators in the context of general ordinary differential equations, whereas in sections 4 and 5 we apply this formalism for constructing the Magnus and Floquet-Magnus expansions, respectively, in the general nonlinear setting. There, we also show how they reproduce the classical expansions for linear differential equations. As a result, both expansions can be considered as the output of appropriately continuous changes of variables rendering the original system into a simpler form. Finally, in section 6 we illustrate the techniques developed here by considering two examples: the Van der Pol oscillator and the nonlinear Schrödinger equation with periodic boundary conditions.

The Magnus expansion for linear systems
There are many ways to get the terms of the Magnus series (3). If we introduce a (dummy) parameter ε in eq. (1), i.e., we replace A by εA, then the successive terms in can be determined by inserting Ω(t) into eq. (1) and computing the derivative of the matrix exponential, thus arriving at [24] εA = d exp Ω (Ω ) ≡ ∞ k=0 1 (k + 1)! ad k Ω (Ω ) = Ω + so that, in terms of equation (8) can be written as Here we have used the notation Now the successive terms R j (t) can be determined by substitution of (10) into (11) and comparing like powers of ε. In particular, this gives Of course, equation (8) can be inverted, thus resulting in , where the B j are the Bernoulli numbers, that is In terms of R, equation (12) can be written as Substituting (10) in eq. (13) and working out the resulting expression, one arrives to the following recursive procedure allowing to determine the successive terms R j (t): Notice in particular that At this point it is worth remarking that any of the above procedures can be used to write each R j in terms of the binary operation £ and the original time-dependent linear operator A, which gives in general one term per binary tree, as in [25,24], or equivalently, one term per planar rooted tree. However, the binary operator £ satisfies, as a consequence of the Jacobi identity of the Lie bracket of vector fields and the integration by parts formula, the so-called pre-Lie relation As shown in [22], this relation can be used to rewrite each R j as a sum of fewer terms, the number of terms being less than or equal to the number of rooted trees with j vertices. For instance, the formula for R 4 can be written in the simplified form If, on the other hand, one is more interested in getting an explicit expression for Ω j (t), ithe usual starting point is to express the solution of (1) as the series where and then compute formally the logarithm of (16). Then one gets [4,30,32,2] with Here σ ∈ S n denotes a permutation of {1, 2, . . . , n}. An expression in terms only of independent commutators can be obtained by using the class of bases proposed by Dragt & Forest [18] for the Lie algebra generated by the operators A(t 1 ), . . . A(t n ), thus resulting in [3] Ω n (t) = 1 n where now σ is a permutation of {1, 2, . . . , n − 1} and d σ corresponds to the number descents of σ. We recall that σ has a descent in i if σ(i) > σ(i + 1), i = 1, . . . , n − 2.

Continuous changes of variables
Our purpose in the sequel is to generalize the previous expansion to general nonlinear differential equations. It turns out that a suitable tool for that purpose is the use of continuous variable transformations generated by operators [17,9]. We therefore summarize next its main features. Given a generic ODE system of the form the idea is to apply some near-to-identity change of variables x −→ X that transforms the original system (20) into where the vector field F (X, t) adopts some desirable form. In order to do that in a convenient way, we apply a one-parameter family of time-dependent transformations of the form such that Ψ 0 (X, t) ≡ X, and x = Ψ 1 (X, t) is the change of variables that we seek.
In this way, one continuously varies s from s = 0 to s = 1 to move from the trivial change of variables x = X to x = Ψ 1 (X, t), so that for each solution X(t) of (21), the function z(t, s) defined by z(t, s) = Ψ s (X(t), t) satisfies a differential equation In particular, we will have that F ( Next, the near-to-identity family of maps X −→ z = Ψ s (X, t) is defined in terms of a differential equation in the independent variable s, by requiring that z(t, s) = Ψ s (z(t, 0), t) for any solution z(t, s) of (23). The map Ψ s (·, t) will be near-to-identity if W (z, t, s) is of the form for some small parameter ε.

Proposition 1 ([17])
Given F and W = εW 1 + ε 2 W 2 + · · · , the right-hand side V of the continuously transformed system (22) can be uniquely determined (as a formal series in powers of ε) from V (X, t, 0) = F (X, t) and where W and V refer to the differentials ∂ x W and ∂ x V , respectively.
Proof. By partial differentiation of both sides in (23) with respect to t and partial differentiation of both sides in (22) with respect to s, we conclude that (24) holds for all x = z(s, t) = Ψ s (x 0 , t) with arbitrary x 0 and all (t, s). One can show that the equality (24) holds for arbitrary (x, t, s) by taking into account that, for given t and s, where are uniquely determined by equating like powers of ε in (25).
In the sequel we always assume that the generator W of the change of variables (i) does not depend on s, and The successive terms in the expansion of V (x, t, s) in Proposition 1 can be conveniently computed with the help of a binary operation £ on maps R d+1 −→ R d defined as follows. Given two such maps P and Q, then P £ Q is a new map whose evaluation at (x, t) ∈ R d+1 takes the value Under these conditions, from Proposition 1, we have that with the notation for the Lie bracket. Equation (27), in terms of where we have used the notation V s (x, t) We thus have the following result: Proposition 2 A change of variables x = Ψ 1 (X, t) defined in terms of a continuous change of variables X −→ z = Ψ s (X, t) with generator and W (x, 0) ≡ x, transforms the system of equations (20) into (21), where f and F are related by and R is given by (29).
Proposition 2 deals with changes of variables such that X = Ψ 1 (X, 0) (as a consequence of W (X, 0) ≡ X), so that the initial value problem obtained by supplementing (20) with the initial condition x(0) = x 0 is transformed into (21) supplemented with More generally, one may consider generators W (·, t) within some class C of timedependent smooth vector fields such that the operator ∂ t : C → C is invertible. Next result reduces to Proposition 2, when one considers some class C of generators W (·, t) such that W (x, 0) ≡ 0, so that ∂ t : C → C is invertible, with inverse defined as within some class C of time dependent smooth vector fields with invertible ∂ t : C → C transforms the initial value problem where f , F , and R = ∂ t W are related by (32), and the binary operator £ : C ×C → C is defined as Notice that the operation £ of (36) satisfies the pre-Lie relation (15), and that this proposition applies, in particular, to the class C of smooth (2π)-periodic vector fields in R d with vanishing average. In that case the operator ∂ t is invertible, with inverse given by

Continuous transformations and the Magnus expansion
Consider now an initial value problem of the form where the parameter ε has been introduced for convenience. As stated in the introduction, the solution x(t) of this problem (20) can be approximated at a given t as the solution y(s) at s = 1 of the autonomous initial value problem This is nothing but the first term in the Magnus approximation of x(t). As a matter of fact, the Magnus expansion is a formal series (31)  This can be achieved by applying Proposition 2 with F (X, t) ≡ 0 and f (x, t) = εg(x, t), i.e., solving and determining the generator W as At this point it is worth analyzing how the usual Magnus expansion for linear systems developed in section 2 is reproduced with this formalism. To do that, we introduce operators Ω(t) and B s (t) such that

Now equation (27) reads
where the binary operation £ defined in (26) reproduces (9). Since B 1 (t) = ε A(t) and B 0 = 0, then (38) is precisely (11). The continuous change of variables is then given by reproduces the Magnus expansion in the linear case. In consequence, the expression for each term W j (x, t) in the Magnus series for the ODE (37) can be obtained from the corresponding formula for the linear case with the binary operation (26) and all results collected in section 2 are still valid in the general setting by replacing the commutator by the Lie bracket (28).

Continuous transformations and the Floquet-Magnus expansion
The procedure of section 3 can be extended when there is a periodic time dependence in the differential equation. In that case one gets a generalized Floquet-Magnus expansion with agrees with the usual one when the problem is linear. As usual, the starting point is the initial value problem where now g(x, t) depends periodically on t, with period T . As before, we apply a change of variables x = Ψ 1 (X, t), defined in terms of a continuous transformation X −→ z = Ψ s (X, t) with generator W = W (x, t) that removes the time dependence, i.e., that transforms (40) into In addition, the generator W is chosen to be independent of s and periodic in t with the same period T . This can be achieved by considering Proposition 2 with F (X, t) := ε G(X; ε) and f (x, t) := ε g(x, t), solving (32) for the series Thus, for the first terms one gets and, in general, where U j only contains terms involving g or the vector fields U m and G m of a lower order, i.e., with m < j. This allows one to solve (42) recursively by taking the average of U j over one period T , i.e., thus ensuring that R j is periodic. Finally, once G and R are determined, W is obtained from (39), which in turn determines the change of variables. If we limit ourselves to the linear case, g( which, with the additional constraint B 1 (t) = A(t), leads to and so on, i.e., exactly the same expressions obtained in [11]. The transformation is now thus reproducing the Floquet-Magnus expansion in the periodic case [11]. Several remarks are in order at this point: 1. This procedure has close similarities with several averaging techniques. As a matter of fact, in the quasi-periodic case, it is equivalent to the high order quasistroboscopic averaging treatment carried out in [15].

2.
A different style of high order averaging (that can be more convenient in some practical applications) can be performed by dropping the condition that W (x, 0) ≡ x, and requiring instead, for instance, that W (x, t) has vanishing average in t.
In that case, the initial condition in (41) must be replaced by X(0) = Ψ −1 1 (x 0 ). The generator W (x, t) of the change of variables and the averaged vector fields ε G(x, t) can be similarly computed by considering Proposition 3 with the class C of smooth quasi-periodic vector fields (on R d or on some smooth manifold) with vanishing average.

Examples of application
The nonlinear Magnus expansion has been applied in the context of control theory, namely in non-holonomic motion planning. The considered systems can be described by equations where dim q > dim u, g i are vector fields and the u i are the controls. The (nonlinear) Magnus expansion allows one to express, locally around a given point in the state space, admissible directions of motions in terms of control parameters, so that the motion in a desired direction can be reformulated as the optimization of those control parameters that steer the non-holonomic system into the desired direction [32,19,21]. In this sense, expression (19) and the corresponding one for the generic, nonlinear case, could be very useful in applications, since it only contains independent terms [20,21].
As mentioned previously, the general Floquet-Magnus can also be applied in averaging. A large class of problems where averaging techniques are successfully applied is made of autonomous systems of the forṁ where h is a smooth function from R n to itself (or more generally from a functional Banach space E space to itself) and A is a skew-symmetric matrix M(R n ) (or more generally a linear operator on E) whose spectrum is included in 2iπ T Z. Denoting x = e −tA u and differentiating leads tȯ so that x now satisfies an equation of the very form (40) with The T -periodicity of g with respect to time stems from the fact that Sp(A) ⊂ 2iπ T Z. For this specific case, relation (42) leads to the following expressions [g(X, σ), g(X, τ ))] dσ, g(X, t) dτ .
If g(x, t) is a Hamiltonian vector field with Hamiltonian H(x, t), then all G i 's are Hamiltonian with Hamiltonian H i 's. These Hamiltonians can be computed through the same formulas with Poisson brackets in lieu of Lie brackets (see e.g. [5]).

Dynamics of the Van der Pol system
As a first and elementary application of previous results, we consider the Van der Pol oscillator, which may be looked at as a perturbation of the simple harmonic oscillator: Clearly, the system is of the form (43) with u = (q, p) T , and is thus amenable to the rewriting (40), where g(x, t) := e −tA h(e tA x) and e tA = cos(t) sin(t) − sin(t) cos(t) .
In short, we haveẋ where V t = − sin(t) cos(t) and ξ t (x) = 1−(cos(t)x 1 +sin(t)x 2 ) 2 (− sin(t)x 1 +cos(t)x 2 ).  Previous formulas give for the first term in the procedure and it is then easy to determine an approximate equation of the limit cycle, i.e. X 2 = 2. As for the dynamics of the first-order averaged system, it is essentially governed by the scalar differential equation on N (X) := X 2 which has two equilibria, namely X 2 = 0 and X 2 = 2. The first one is unstable while the second one is stable. However, the graphical representation (see Figure 1) of the solution of (47) soon convinces us that the true limit cycle is not a perfect circle.
In order to determine a better approximation of this limit cycle, we thus compute the next term of the average equation from formula (45): where D(X) = 32 − 24 X 2 2 + 5 X 2 4 − 88 X 1 2 + 21 X 1 4 + 10 X 1 2 X 2 Now, considering the new quantity L(X) = N (X) + εQ(X) with Q(X) = νX 1 X 3 2 , we see from the second-order averaged equatioṅ for ν = − 1 2 . A more accurate description of the limit cycle is thus given by the equation

The first two terms of the averaged nonlinear Schrödinger equation (NLS) on the d-dimensional torus
Next, we apply the results of Section 5 to the nonlinear Schrödinger equation (for a introduction to the NLS see for instance [8,14,13]) is our working space. We hereafter conform to the hypotheses of [12] and assume that h is a real-analytic function and that s > d/2, ensuring that E is an algebra 1 . The operator −∆ is self-adjoint non-negative and its spectrum is so that by Stone's theorem, the group-operator exp(it∆) is periodic with period T = a 2 2π . We may thus rewrite Schrdinger equation as we did with equation (43). However, we shall instead decompose ψ(t, z) = q(t, z) + ip(t, z) in its real and imaginary parts and derive the corresponding canonical system in the new unknown that is to saẏ where we have denotedu = ∂ t u, u 2 R 2 = (u 1 ) 2 + (u 2 ) 2 = (u, u) R 2 (u 1 and u 2 are the two components of u), J is given by (48) and The operator D if self-adjoint on L 2 (T d a )×L 2 (T d a ) and an obvious computation shows e tJ −1 D = cos(t∆) sin(t∆) − sin(t∆) cos(t∆) so that e tJ −1 AD is a group of isometries on L 2 (T d a ) × L 2 (T d a ) as well as on H s (T d a ) × H s (T d a ). Owing to (49) it is furthermore periodic (for all t, R t+T = R t ), with period T = a 2 2π . The very same manipulation as for the prototypical system (43) then leads tȯ Now, it can be verified that , where x 1 1 and x 2 1 are the two components of x 1 and similarly for x 2 . Hence, by definition of the gradient, we have that Furthermore, and ∀(x 1 , x 2 , x 2 ) ∈ E 3 (J∂ x g(x 1 , t)x 2 , x 3 ) = (x 2 , J∂ x g(x 1 , t)x 3 ).
Finally, if φ 1 and φ 2 are hamiltonian vector fields, with hamiltonians Φ 1 and Φ 2 , then where the Poisson bracket is defined by Now, the first term of the averaged vector field G(X, ε) is simply In order to obtain the second term, we use the simple fact that for any δ ∈ H s (T d a ) the derivatives w.r.t. x in the direction δ may be computed as Inserted in the expression of G 2 we thus obtain the following expression for the ε 2term of the averaged equation [g(X, τ ), g(X, t)] dτ = I 1 (X) + I 2 (X) with As already mentioned, both G 1 and G 2 are Hamiltonian with Hamiltonian H 1 and H 2 which could have been equivalently computed from H(x, t) in (51) (see Remark 1).