High order linearly implicit methods for evolution equations

This paper introduces a new class of numerical methods for the time integration of evolution equations set as Cauchy problems of ODEs or PDEs. The systematic design of these methods mixes the Runge-Kutta collocation formalism with collocation techniques, in such a way that the methods are linearly implicit and have high order. The fact that these methods are implicit allows to avoid CFL conditions when the large systems to integrate come from the space discretization of evolution PDEs. Moreover, these methods are expected to be efficient since they only require to solve one linear system of equations at each time step, and efficient techniques from the literature can be used to do so. After the introduction of the methods, we set suitable definitions of consistency and stability for these methods. This allows for a proof that arbitrarily high order linearly implicit methods exist and converge when applied to ODEs. Eventually, we perform numerical experiments on ODEs and PDEs that illustrate our theoretical results for ODEs, and compare our methods with standard methods for several evolution PDEs.


Introduction
The goal of this paper is to introduce a new class of methods for the time integration of evolution problems, set as (systems of) deterministic ODEs or PDEs.This class consists in methods of arbitrarily high order that require only the solution of one linear problem at each time step: no nonlinear system is to be solved.As is usual in the litterature, we call these methods linearly implicit.They rely on the combination of a classical collocation Runge-Kutta method with a specific treatment for the nonlinearity.
In particular, we show that, using the methods developed in this paper, one can solve numerically virtually any ODE, up to any order, by solving only linear systems at each time step.Moreover, we believe that this new class of methods can help dramatically reducing the computational time in several cases of time integration of evolution PDEs.Indeed, after space discretization of an evolution PDE, if one uses, say, an implicit (for stability reasons) Runge-Kutta method, then one needs to solve a nonlinear problem in high dimension at each time step (and one may use a fixed-point method or a Newton method to do so, for example).With schemes for dissipative problems with gradient flow structure.The auxiliary variable in this context is used to ensure discrete energy decay.The order of the methods (1 or 2 in the references above) is not the main issue.
Let us mention two additional goals that the authors aim at tackling with the methods introduced in this paper.First, the authors would like to be able to develop a stability (and convergence) analysis for stiff problems, i.e. an analysis with constants that depend only on the class of the linear part of the vector field (later referred to as , see (2.1)) and not on that linear part itself.This would allow for the numerical treatment of evolution PDE problems, as well as their space discretizations.This will be achieved in a forthcoming work, even if numerical examples of PDE problems are presented in Section 3. Second, the authors would like to build up high order linearly implicit methods with suitable qualitative properties (e.g.energy decay for dissipative problems, or energy preservation for hamiltonian problems).For example, the relaxation method [5], which belongs to the class of linearly implicit methods described in this paper, preserves an energy when applied to the nonlinear Schrödinger equation.For this reason, this paper only deals with constant time step methods.
The outline of this paper is as follows.In Section 2, we introduce the methods for a general semilinear evolution problem (ODE or PDE) and we introduce specific notions of stability, consistency and convergence for our class of methods.Morover, we show, in a constructive way, that stable methods of arbitrarily high order exist in Theorem 2.5 and Corollary 2.6.The main theoretical result of this paper is that one can build up arbitrarily high order convergent linearly implicit methods for ODEs (Thm.2.9).We conclude Section 2 with examples of methods of order 1, 2, 4 and 6.In Section 3, we provide numerical examples of solutions of ODEs and PDEs.These numerical experiments illustrate the convergence result of Theorem 2.9 for evolution ODEs.Moreover, they indicate that the result of Theorem 2.9 is still valid in several PDE contexts.We consider for example a NLS equation in 1d and 2d and nonlinear heat equation in 1d.The main result of the numerical experiments of Section 3 is that, for ODEs, the linearly implicit methods do not dramatically outperform classical methods from the literature with the same order, no matter whether they are implicit or explicit (see Sect. 3.1).However, for the approximation of evolution PDEs in 1d (see Sect. 3.2.1),with moderate space discretization, the linearly implicit methods show performances comparable to that of explicit methods.Moreover, for the approximation of evolution PDEs in 2d (see Sect. 3.2.2) with precise space discretization (leading to high number of unknowns), the linearly implicit methods developed in this paper manage to outperform standard methods from the literature with the same order.

Introduction of the methods
We consider a semilinear autonomous evolution equation of the form where  is a linear differential operator and  is a nonlinear function of .One can think for examples of the NLS equation, the nonlinear heat equation, or a simple ODE (see Sect. 3 for actual examples).We start at time  " 0 with an initial datum  0 (with 0 in superscript) in some functional space so that the Cauchy problem is well-posed on some interval r0,  ‹ q with  ‹ ą 0. We choose ℎ ą 0 and set   " ℎ for  P N as long as   ă  ‹ .
Let us now start with the presentation of the new class of methods.Assume a collocation Runge-Kutta method with  ě 1 stages is given with coefficients 0 ď  1 ă ¨¨¨ă   ď 1, p , q 1ď,ď and p  q 1ďď .We denote by  the vector p  q 1ďď and by 1 the vector of size  with all entries equal to one.
We denote by  the exact solution of (2.1) and we set pq "  pp, ¨qq.We assume we are given  approximations  ´1` " p ´1 ` ℎq 1 ď  ď , and another approximation   " p  , ¨q.For possible ways of computing the  first approximations, we refer to Remark 2.11.For the numerical initial datum  0 (with subscript 0), we consider an approximation of the exact initial datum  0 (with superscript 0) of equation (2.1) or its exact value.For p 1 , . . .,   q P R  and  P ℳ  pRq to be chosen later, we define explicitly p `1 , . . .,  ` q with the relation Then, we define, linearly implicitly p ,1 , . . .,  , q as the solution of the Runge-Kutta like system Last, we set explicitly p `` q , . (2.4) The steps (2.2)-(2.4)define a linearly implicit method p `1 ,  `1 , . . .,  ` q " Φ ℎ p  ,  ´1`1 , . . .,  ´1` q.
In order to achieve order 2 with the relaxation method, C. Besse introduced a single auxiliary unknown  on a staggered grid, corresponding to the relation  " || 2 .Since we want to achieve higher orders, we decide to introduce several auxiliary unknowns on a staggered grid with  points, corresponding to the relation  "  pq.
Note that the convergence of order 2 of the relaxation method for the NLS equation is a difficult result and is not a consequence of the results of this paper, which only deals with ODEs in the theoretical part (Sect.2) and allows for PDEs for illustration purposes (Sect.3).Indeed, proving the convergence of a numerical time integration method applied to a PDE requires a functional analysis framework adapted to the PDE at hand and cannot in general be done once and for all.For the relaxation method applied to the NLS equation, the convergence of order 2 is proved in [7].

Consistency and stability of the step (2.2)
Let us denote by pq the spectral radius of the matrix , i.e. the biggest modulus of its complex eigenvalues.In view of relation (2.2), we decide to set the following definitions for the stability and consistency of the step (2.2).
for some norm on ℳ  pRq.The step (2.2) is said to be strongly stable if pq ă 1.
Remark 2.3.In the definition of the stability above, the boundedness of the sequence p  q ě0 is independant of the norm chosen on ℳ  pRq.Moreover, it is equivalent to the fact that pq ă 1 or pq " 1 with simple Jordan blocks for  for all eigenvalues of modulus 1.In particular, if the step (2.2) is strongly stable, then it is stable.The converse is not true in general.For example, the classical relaxation method of Remark 2.1 is stable but not strongly stable.
In order to define the consistency of step (2.2), we introduce the  ˆ square matrices   ,  ´1 and Θ defined by -for all  ě 1 and  ě 1, ` ℎ  ˘ " p  ℎq ´1 , -for all  ě 1 and  ě 1, ` ℎ ´1 ˘ " pp  ´1qℎq Since Θ is zero except maybe on its first column, we have Θppℎqq ´1 " Θ.The definition of the step (2.2) of the method (2.2)-(2.4)depends on the  2 coefficients of the matrix  and the  coefficients  1 , . . .,   .Requiring that the step (2.2) is of order  provides us with  2 linear equations between these unknowns (see relation (2.7)).In Theorem 2.5, we prove that we can add  equations involving these unknows by imposing the spectrum of the matrix  and that the system that we obtain has indeed a unique solution.This will allow in particular to prove the existence of stable and strongly stable steps (2.2) with order  (see Cor. 2.6).
Theorem 2.5.Assume  1 , . . .,   are fixed and distinct as above.For all disctinct  1 , . . .,   P Czt1u, there exists a unique pp 1 , . . .,   q, q P C  ˆℳ pCq such that the step (2.2) is of order  and the spectrum of the matrix  is exactly t 1 , . . .,   u.If moreover the set t 1 , . . .,   u is stable under complex conjugation, then p  q 1ďď P R  and  P ℳ  pRq.
Proof.We set  " ` 1 ´1 ˘´1  1   .Note that  is in fact independant of the choice of the p  q 1ďď .Indeed, it is the matrix of the linear mapping  pq Þ Ñ  p `1q in the canonical basis of R ´1 rs.This means that the coefficients p  q 1ď,ď of  are given by   " 0 if  ă  and   " `´1 ´1 ˘otherwise.In particular, it is upper triangular and its diagonal elements are equal to 1. Assuming step (2.2) is consistent of order , with (2.7), we obtain so that the matrix  is similar to  ´ , where  is the matrix ` 1

𝑐´1
˘´1 Θ.Note that all the coefficients of  are equal to 0, except maybe on the first column.We shall denote by  1 , . . .,   the coefficients in the first column of  and by  1 the first column of  .Given distinct  1 , . . .,   P Czt1u, the existence and uniqueness of  and Θ such that step (2.2) has order  and the spectrum of  is exactly t 1 , . . .,   u is equivalent to the existence and uniqueness of  1 , . . .,   P C such that  ´ has spectrum t 1 , . . .,   u.
Let us fix  P t1, . . ., u.The existence of an eigenvector for  ´ for the eigenvalue   is exactly the existence of a nontrivial vector   P C  such that p ´ q  "     .Let us denote by  the identity matrix of size  and  the upper triangular matrix such that  "  ´ .Observe that, in view of the definition of  , the entries above the diagonal of  are negative.‰ 0 and we can impose without loss of generality that  pq 1 " 1. Together with (2.9), this implies, by recursion, that (2.10) The equation on the first line in (2.10) is of the form where for all  P t1, . . ., u,   is a polynomial of degree exactly  (remind that the matrix  is strictly upper triangular with negative entries above the diagonal).Moreover, if the relation (2.11) is verified for some p 1 , . . .,   q, then there exists a solution   of (2.9) with  pq 1 " 1: one just has to compute the components of   in (2.10) one after the other to obtain an eigenvector of  ´ for the eigenvalue   .As a summary, we have proved that, for all  P t1, . . ., u,   is an eigenvalue of  ´ if and only if (2.11) holds.Therefore, the fact that the spectrum of  ´ is t 1 , ¨¨¨,   u is equivalent to the linear system » ---- (2.12) Since for all  P t1, . . ., u, the polynomial   has degree  and the p  q 1ďď P pCzt1uq  are distinct, the system above is invertible, so that it has a unique solution.Since (2.12) has a unique solution in C  and the polynomials p  q 1ďď have real coefficients, it is easy to check that, if the set t 1 , . . .,   u is moreover stable under complex conjugation, then  1 , . . .,   are real numbers, and so are  1 , . . .,   and the matrix  has real coefficients using (2.8).This proves the theorem.
-There exists  P ℳ  pRq and  1 , . . .,   P R such that the step (2.2) is strongly stable and has order .
Proof.Choose  1 , . . .,   P Czt1u distinct such that for all , |  | ď 1 to obtain a stable method (or for all , |  | ă 1 to obtain a strongly stable method) in such a way that the set t 1 , . . .,   u is stable under complex conjugation and apply Theorem 2.5.
Remark 2.7.In order to actually build the matrix  and the coefficients  1 , . . .,   that define step (2.2) so that this step is stable (respectively strongly stable) and has order , it is sufficient to fix distinct  1 , . . .,   as above, and choose distinct  1 , . . .,   P Czt1u with modulus less (respectively strictly less) than 1, and in such a way that the set t 1 , . . .,   u is stable under complex conjugation.Then, one forms system (2.12), the rows of which are the first rows of the right-hand side of (2.10) for different values of   , and one solves it for  1 , . . .,   .One easily computes Θ from  using the fact that Θ "  1 ´1  .In the end, the matrix  is computed using (2.7).Examples are provided in Section 2.4.

Convergence of the method (2.2)-(2.4)
In this section, we prove that the methods presented above, provided that they involve a step (2.2) with strong stability and order , and a Runge-Kutta collocation method of order at least , are indeed convergent in finite time, with order , when applied to an ODE with sufficiently smooth vector field.We assume the unknown  of equation (2.1) is scalar and that  " 0. In fact, up to a change of unknown, any ODE with an equilibrium can be cast into this form: 1 pq "  ppqqpq. (2.13) Our methods and results extend to systems of ODEs of the form (2.13) where the unknown  is vector-valued and  is a given smooth matrix-valued function.Similarly, our methods and results extend to the case of complex-valued functions.But for the sake of simplicity we focus on the real valued scalar case.We assume  is defined and smooth on some open subset Ω of R. We fix  0 P Ω.There exists a unique maximal solution to the Cauchy problem (2.13) for p0q "  0 .This solution is defined on an open interval of the form p ‹ ,  ‹ q with ´8 ď  ‹ ă 0 ă  ‹ ď `8.We fix  P p0,  ‹ q.Since the maximal solution is smooth, we have sup Pr0, s | ppqq| ă `8.Since it is defined on the compact interval r0,  s, one can choose  ą 0 such that  " tpq ` |  P r0,  s, We set  " sup Pr0, s | ppqq| and we choose  ą 0 such that  ` ą sup P | pq|.
We discretize the time as in Section 2.1, with ℎ small enough to ensure that  ‹ ă  ´1.We start by focusing on the consistency of the method.Namely, we set for all  P N such that   ď  , Similarly, we define  2  as the vector of R  with entry number  equal to and (2.17) Lemma 2.8.Assume that the function  is sufficiently smooth,  0 P Ω and  P p0,  ‹ q.Suppose moreover that the numerical coefficients p  q 1ďď , p , q 1ď,ď and p  q 1ďď define a Runge-Kutta collocation method of order  and that the step (2.2) is of order .For any norm on R  , there exists a constant  ą 0 such that for a sufficiently small ℎ ą 0, max ě0, `1ď Since the step (2.2) is of order , we have using (2.6) Moreover with (2.7) we have for some  ą 0 and all sufficiently small ℎ ą 0. This proves (2.18).Since the Runge-Kutta method with coefficients  , and   is a collocation method of order at least  at points   , the bounds (2.19) and (2.20) are classical (see e.g.[20], Sect.II.1.2).
Theorem 2.9.Assume that the function  is sufficiently smooth,  0 P Ω and  P p0,  ‹ q.Suppose moreover that the numerical coefficients p  q 1ďď , p , q 1ď,ď and p  q 1ďď define a Runge-Kutta collocation method of order  and that the step (2.2) is strongly stable and of order .Provided that  is fixed as in (2.14) and  and  accordingly, there exists constants  ą 0 and ℎ 0 ą 0, such that, for all ℎ P p0, ℎ 0 q, if the initial data  0 P Ω (respectively p ´1`1 , ¨¨¨,  ´1` q P  p q  ) is sufficiently close to its exact analogues  0 P Ω (respectively p pp ´1 `1 ℎqq, ¨¨¨,  pp ´1 ` ℎqqq P  p q  ) in the sense of relations (2.38) and (2.39), then for all  P N such that  ´1 ď  , the step (2.3) has a unique solution in R  and for all  P N such that   ď  , Let us first introduce all the notations we use in the proof.We denote by Γ  the vector of R  with component  equal to  ` .Let us define the convergence errors   P R  with component number  equal to p  q  "  pp  ` ℎqq ´` ,   P R  with component number  equal to  , " p  ` ℎq ´, (provided  , is well defined), and   P R with   " p  q ´ .We set   " max 0ďď |  |.We denote by | ¨|8 the norm on R  defined as the maximum of the absolute values of the components of the vectors.Moreover, we denote by } ¨}8 the norm on ℳ  pRq induced by | ¨|8 .In the following proof, the letter  denotes a positive real number which does not depend on ℎ (but depends on  and  in particular) and whose value may vary from one line to the other.
Proof.Since step (2.2) is strongly stable we have pq ă 1.Therefore, there exists a norm | ¨| on R  such that the norm } ¨} induced by this norm on ℳ  pRq satisfies }}  ă 1.In the following we set  " }}  .Since R  is of dimension , there exists a  P p0, 1s such that for all  in R  , ||  ď || 8 ď 1  ||  .We divide the proof in two parts.First we assume an a priori bound for the numerical solution.Namely we assume that for all  such that   ď  : 3) has a unique solution p , q 1ďď in R  , -(H3)   P  .
We show that, in this case, we have an explicit bound for the convergence errors   and   (see Eqs. (2.35) and (2.37)).
Second, we assume that ℎ 0 and the initial errors  ´1 and  0 are small enough and we show that the bounds of the first part of the proof are indeed satisfied.
First part.In addition to the bounds above, we assume that ℎ P p0, 1q and  satisfy  `1 ď  .Substracting (2.2) from (2.15) we obtain We infer that where the constant  is the product of the Lipschitz constant of  over the compact  times |p Using the norm equivalence and (2.30) in (2.29) we obtain which gives with Lemma 2.8 Using the maximal error defined previously and the fact that  ă 1 since the step (2.2) is strongly stable, we have and then Therefore by triangle inequality we have And then, the hypothesis (H1) of the first part is satisfied for  " 0.Moreover, with (2.38), we have | 0 ´0 | ď  so that  0 P  and the hypothesis (H3) of the first part is satisfied with  " 0. We infer that the system (2.3) (with  " 0) has a unique solution in R  since we assumed ℎ ď ℎ 0 ă 1{p}} 8 p `qq.This implies that the hypothesis (H2) of the first part is satisfied for  " 0. Then we can apply the analysis of the first part to obtain (2.35) and (2.37) with  " 1.Using (2.38) and (2.39), we infer that hypotheses (H1)-(H3) are satisfied with  " 1 and the result follows by induction on .
Remark 2.10.The proof may look like following the usual strategy for the proof of convergence of a numerical method.However, note that the definition of strong stability (Def.2.2) plays a central role in the proof, and this is not the case in classical theory.Moreover, estimations such as (2.30) in (2.29) to obtain (2.31) are not classical.
Remark 2.11.Given  0 P Ω close to  0 P Ω, before starting the time-stepping method (2.3) and (2.4), one needs to first compute the  approximations p ´1` q 1ďď of p ppp´1 ` qℎqqq 1ďď .These quantities enter the error estimate (2.25).This computation can be done efficiently by determinating  approximations of ppp´1 ` qℎqq 1ďď using a sufficiently high order method for (2.1) backwards in time, provided the equation makes sense.Alternatively, for example for the nonlinear heat equation, for which running (2.1) backwards in time makes no sense, one can use a sufficiently high order method from  0 to compute  `1 approximations of pp  ℎqq 1ďď and pℎq over one time step, and then start the linearly implicit method (2.3) and (2.4) from time ℎ until time  .

Examples of linearly implicit methods
In this section we present possible choices of methods of order 1, 2, 4 and 6.The general building procedure is the following: We choose  P N ‹ , we fix 0 ď  1 ă  2 ă ¨¨¨ă   ď 1 and we compute  , and   for 1 ď ,  ď  using the formulas where ℒ  p q " ś  "1 ‰ p ´ q p  ´ q is the th Lagrange polynomial at points  1 , . . .,   .This way, the coefficients  , ,   and   are those of a Runge-Kutta collocation method.Next we choose  1 , . . .,   P Czt1u with moduli strictly less than 1, all distincts and in such a way that the set t 1 , . . .,   u is invariant under complex conjugation.We compute the polynomials  1 , . . .,   appearing in (2.11) defined using (2.10) in the proof of Theorem 2.5.We solve (2.12) for  1 , . . .,   and compute  1 , . . .,   using Θ " ` 1 ´1 ˘ .Finally, we compute the matrix  using (2.7).This way, we define a step (2.2) that is strongly stable and of order  (see Defs. 2.2 and 2.4).Using Theorem 2.9, the numerical method (2.2)-(2.4) is convergent of order .
Choosing   "  p´1q  3 2 for  " 1, . . ., 6, we have Remark 2.12.In the linearly implicit methods given as examples above, the choice of p  q 1ďď is somehow arbitrary, and the only condition we impose is that they ensure that step (2.2) is strongly stable (see Def. 2.2).This implies that these methods are convergent for ODEs (see Thm. 2.9).In order to ensure additional features of a linearly implicit method, this choice has to be made carefully (see e.g.Sect.3.3 a linearly implicit method that preserves non-negativity and energy-dissipation for a nonlinear heat equation).For a general class of evolution PDE, the choice of the collocation method as well as that of p  q 1ďď has to be made carefully, in particular to ensure a taylored stability property of the linearly implicit method.This question will be investigated in a forthcoming paper.

Numerical experiments
In this section, we illustrate the properties of the methods described in Section 2.1 and analysed in Section 2.3.We first present numerical examples on ODEs, with a scalar case in Section 3.1.In particular, we illustrate the results above, such as Theorem 2.9, for several methods introduced above, and we compare the results we obtain with that obtained using other classical numerical methods.Then, we present numerical experiments for PDEs that fit the framework used in Section 2.1 but do not fit stricto sensu the framework of the analysis carried out in Section 2.3.This allows for comparison with classical methods for the same problems anyway.We first focus on a nonlinear Schrödinger equation in Section 3.2 and then move to a nonlinear heat equation in Section 3.3.The methods described and analysed in this paper would also be relevant for several other examples of semilinear evolution equation of the form (2.1).
When comparing the efficiency of numerical methods in this section, we consider as a measure of performance the (lowest possible) CPU time required to achieve a given precision on the numerical result.This CPU time has indeed disadvantages since it depends on the algorithms used to solve the problems, the software used to implement the algorithms and the machine on which the software is run.However, we believe one cannot talk about efficiency whithout taking into account some form of CPU time.And, for reproducibility issues, we detail below as much as possible which discretizations and algorithms are used to implement the numerical methods that we consider.Moreover, we try to be as fair as possible when implementing methods from the litterature to compare them with the linearly implicit methods introduced in this paper.
As we shall see in this section, the efficiency of the linearly implicit methods introduced in this paper is similar to that of classical one-step methods with constant step size from the literature (see Sect. 3.1).
In contrast, the linealy implicit methods sometimes outperform standard methods when applied to several evolution PDE problems, that we consider, once discretized, as high dimensional systems of ODEs (see Sects.3.2 and 3.3).In the following, the computations are carried out using MATLAB and the linear systems are solved using the backslash MATLAB command.In particular, we do not build a taylored method to solve the linear systems numerically and the gain in computational time one can obtain using linearly implicit methods can surely be improved using taylored methods depending on the matrix structures.This choice is not optimal in terms of efficiency, but is fairly similarly done for all the methods below.

Application to a scalar nonlinear ODE
We consider the scalar ODE  1 pq " ´pq ´2 pq.(3.1)This corresponds to taking  as minus the identity operator and  pq " ´ in (2.1).The exact maximal solution starting from  0 ą 0 at  " 0 is given for  ě 0 by Ẅe start with methods of order 1.We use the linearly implicit method of order 1 introduced in Section 2.4.We compare the results we obtain on the problem above with the Euler implicit and explicit schemes as well as the Lie splitting method.We choose  0 "  0 " 1{3,  ´1`1 "  0 "  p 0 q and the final time  " 2. The results are displayed in Figure 1.The global error is defined as   (with the notations of the proof of Thm.2.9) at final time  with  such that  ℎ "  , where ℎ is the time step.Numerical experiments indicate that the four schemes are of order 1.For the linearly implicit scheme, this is a consequence of Theorem 2.9.Moreover the CPU time required to reach a given numerical error is much lower for the Lie splitting than for the linearly implicit method and for the linearly implicit method than for the explicit Euler scheme and the implicit Euler scheme.We then consider methods of order 2. We compare the linearly implicit method of order 2 defined in Section 2.4 for Gauss points with other methods of the literature: the midpoint method with Butcher tableau 1{2 1{2 1 , the RK2 method with Butcher tableau 0 0 0 , and the Strang splitting method.We choose  0 "  0 " 0.9,  ´1`1 "  ppp´1 `1 qℎqq,  ´1`2 "  ppp´1 c2 qℎqq and the final time  " 2. The results are displayed in Figure 2. Once again the four methods are of order 2.This is a consequence of Theorem 2.9 for the linearly implicit method.The CPU time required for a given numerical error is much lower for the Strang splitting scheme than for the other three methods which perform similarly.
Similar results are obtained (but not displayed here) for methods of order four and six which illustrate Theorem 2.9 for the corresponding linearly implicit methods introduced in Section 2.4.Moreover the CPU time required to reach a given numerical error is always higher for the linearly implicit schemes than for other classical methods of the same order for the ODE (3.1).

One dimensional nonlinear Schrödinger equation: The soliton case
In this section, we consider the nonlinear one dimensional Schrödinger equation: which corresponds to the evolution problem (2.1) with  " B 2  and  pq " || 2 .We consider the initial condition where  ą 0 and  "  2 {16, so that the corresponding exact solution of (3.2) is the zero speed soliton and reads p, q " We use  " 4 and  " 1 for the numerical simulations.The final time is set to  " 5.For the space discretization, we consider the interval r´50, 50s with homogeneous Dirichlet boundary conditions since the exact solution (3.3) decays very fast when || tends to `8.We use 2 14 equispaced points in space for methods of order 1 in time and 2 18 equispaced points in space for methods of order 2 in time.For splitting methods, one has to integrate numerically equation (3.2) with  " 0. This is done via the approximation where  denotes the identity matrix,  the Laplacian with Dirichlet boundary conditions matrix, and ℎ ą 0 the small time step.In particular, the splitting methods we use below are also linearly implicit (the nonlinear part of the equation is integrated exactly).The numerical error we consider for all numerical methods is the discrete  2 -norm of the difference between the numerical solution and the projection of the exact solution (3.3) on the space grid at final time  .First we compare methods of order one.We consider the linearly implicit method of order 1 introduced in Section 2.4, the implicit Euler scheme and the Lie splitting scheme.For the linearly implicit method, we initialize the scheme with  0 "  0 " p0, ¨q,  ´1`1 "  ppp´1 `1 qℎ, ¨qq.The results are given in Figure 3.The figure on the left hand side shows that the three methods are of order 1.This illustrates the fact that the conclusion of Theorem 2.9 for the linearly implicit method holds numerically in this PDE context.For such a simulation, we can see, on the figure on the right hand side, that now the CPU time required to reach a given error is smaller for the linearly implicit method than for the fully implicit Euler method.However the Lie splitting method is the least time consuming method since it is explicit (in fact our implementation of the Lie splitting method makes it linearly implicit, see (3.4)) and has a good error constant.
We then consider methods of order 2. For this experiment we use the linearly implicit method of order 2 introduced in Section 2.4 for the uniform points, the Crank-Nicolson scheme [13,15] (for which we solve the nonlinear system using a fixed point algorithm) and the Strang splitting method [27].For the implementation of the Strang splitting method, we use the classical conjugation with the Lie splitting method which we recall briefly below and relies on the identity where Φ ℎ denotes the numerical flow of the nonlinear part of (3.2) defined componentwise using the function  Þ Ñ exppℎ|| 2 q, and  is any nonnegative integer.The numerical flow of the Lie splitting method is Φ Lie ℎ " Φ ℎ ˝exppℎq and that of the Strang splitting method is Φ Strang ℎ " exppℎ{2q ˝Φℎ ˝exppℎ{2q.Therefore, relation (3.5) also reads ´ΦStrang for all nonnegative integer .The implementation of the Strang splitting method we use for numerical simulations uses both the right hand side of the equation (3.6) and the approximation formula (3.4).For the linearly implicit method, we initialize the scheme with  0 "  0 " p0, ¨q,  ´1`1 "  ppp´1 `1 qℎ, ¨qq,  ´1`2 "  ppp´1 `2 qℎ, ¨qq.The results are displayed in Figure 4.The numerical order of each method is the one expected i.e. 2. This illustrates once again the numerical relevance of Theorem 2.9 beyond the ODE context.
As we can see on the figure on the right hand side, the CPU time required to reach a given error for the Crank-Nicolson method is higher than the one for the linearly implicit method of order 2 which is also a little higher than the one for the Strang splitting method.
Remark 3.1.For this example, if instead of using the linearly implicit method of order 2 with uniform points, we use the linearly implicit method of order 2 with Gauss points, we numerically obtain the superconvergence of the method and a numerical order equal to 4. This is due to the fact that in this particular case the modulus of the solution is constant in time as it can be seen on (3.3).

Two dimensional nonlinear Schrödinger equation
In this section, we consider the following 2D nonlinear Schrödinger equation:  with homogeneous Dirichlet boundary conditions on the domain represented in gray in the Figure 6 with   "   " 1,   " 2 and   " 3. The initial datum we chose reads  0 p, q " sinp2q sinp2q expp2q, (3.8) when p, q belongs to the domain.In this particular case, the spectrum of the Laplace operator is not accessible and one cannot use efficiently spectral methods such as exponential Runge-Kutta methods or Lawson methods [6].
We use a finite differences discretization in space with, for  P N ˚,    `1 points in the -direction and    `1 points in the -direction in such a way that the step is the same in the two directions.This way, the numerical unknown   at time   is a vector of  " pp  ´1q ´1q ˆp   ´1q ` ˆpp  ´1q ´1q complex numbers.No matter the way we label the unkowns, the matrix  of the Laplace operator with homogeneous  Dirichlet boundary conditions is a sparse matrix of size  ˆ .For the numerical simulations we use  " 50 which gives  " 12 251 unknowns.Moreover, we consider  " 0.5 as a final time.The methods we consider are the two linearly implicit methods of order 2 introduced in Section 2.4 (one with Gauss points, the other one with uniform points), the Crank-Nicolson scheme and the Strang splitting method.As one has no direct access to the exact solution of (3.7) with initial condition (3.8), we precompute as a reference solution the numerical solution provided by a Runge-Kutta method at Gauss points with 5 stages (which therefore has order 10) with a time step of 10 ´3.We initialize the linearly implicit methods with  ´1`1 and  ´1`2 computed using one step of a backward Crank-Nicolson scheme.Our numerical results are displayed in Figure 7.As expected, the order of each method above is 2. For the two linearly implicit methods, this again illustrates that the results of Theorem 2.9 extend numerically to this PDE case.Note that, the constant of order is really better for the linearly implicit method with Gauss points than all the other ones.Moreover the CPU time required to achieve a given precision is also smaller for the linearly implicit method with Gauss points.This is the first example where a linearly implicit method developped in this paper clearly outperforms implicit as well as explicit standard methods from the literature.Remark 3.2 (Preservation of mass in NLS equations by linearly implicit methods).If we consider a linearly implicit method defined by a Runge-Kutta collocation method of order  satisfying the Cooper condition @p, q P t1, ¨¨¨, u 2 ,     "    , `  , , ( that is to say a Runge-Kutta collocation method at Gauss' points, then this linearly implicit method preserves the mass (i.e. the squared  2 -norm) for the NLS equation, regardless of the physical dimension of the problem.This property does not depend on the choice of the matrix  and vector p 1 , ¨¨¨,   q in step (2.2) as long as they are real-valued.This preservation property relies on the fact that, provided the initial p ´1` q 1ďď are purely imaginary, so will be all the p ` q 1ďď .This property will be detailed further in a forthcoming paper dealing with the time integration of PDEs using linearly implicit methods.For example, the linearly implicit methods of order 1, 4 and 6 from Section 2.4 do not satisfy the Cooper condition (3.9), hence they do not preserve the mass for the NLS equation.In contrast, the first method of order 2 of Section 2.4, which uses Gauss' points, preserves the mass of the solution of the NLS equation (and the second method of order 2 does not).

Application to the nonlinear heat equation
In the previous sections, we have proved and illustrated that the linearly implicit methods developed in this paper have good quantitative properties.The goal of this section is to illustrate that they can indeed also have good qualitative properties.Indeed, on some nonlinear heat equation with gradient-flow structure, we give an example below of a linearly implicit fully discrete scheme which preserves the nonnegativity of the solution (just as the exact flow does) as well as the decay of a discrete energy which is consistent with the continuous energy of the problem.
Let us consider the one dimensional nonlinear heat equation given by Finally if  0 ě 0 on Ω then for all  P r0,  ˚q, pq ě 0 on Ω.
Let us denote by  the number of unknowns, so that  " 100{p `1q.We denote by x¨, ¨y the scalar product on R  defined for ,  P R  by x, y "  ř  "1 pqpq and by } ¨}2 the associated norm.Moreover, for all  P R  , we denote by  ˝2 the vector of R  with component  equal to  ˝2pq " pq 2 .
Then, the stage (2.2) reads here and the stages (2.3), and (2.4) can be summarized by where  denotes the matrix of the Laplacian operator with homogeneous Dirichlet boundary conditions on Ω on the equispaced grid, with space step size , as defined after (3.4).We still denote by  0 the evaluation of the initial datum  0 on the equispaced grid.In addition, we choose for  ´1{2 the evaluation of  p 0 q on the same grid.
The fully discrete energy associated to the numerical scheme is defined for ,  P R  by Note that this formula is consistent with the continuous energy  defined in (3.11).
Theorem 3.4.Let us assume  0 P  1 0 pΩq.Still denote by  0 the projection of  0 onto the equispaced grid with  interior points.Choose  P p0,  ˚q.
(1) Let us assume that there exists ℎ 0 ,  0 ą 0 such that for all ℎ P p0, ℎ 0 q and all  P p0,  0 q with ℎ ă  2 , the sequence p `1{2 q ě0 is bounded in R  with the maximum norm as long as p `1{2qℎ ď  .Then, for all ℎ P p0, ℎ 0 q,  P p0,  0 q and  such that p `1{2qℎ ď  and ℎ{ 2 ă 1,  `1{2 is a nonnegative real-valued vector.Moreover assuming  0 ě 0, there exists a constant ℎ 1 P p0, ℎ 0 q such that for all ℎ ď ℎ 1 and  P p0,  0 q with ℎ{ 2 ă 1, the sequence p  q ě0 is a sequence of nonnegative vectors as long as ℎ ď  .(2) For all ℎ P p0, ℎ 0 q and for all  P N such that p `1qℎ ď  , the sequence p  ,  ´1{2 q ě0 satisfies   p `1 ,  `1{2 q ď   p  ,  ´1{2 q. (3.16) Proof.We will use  ,  and  matrices as defined for example in the Chapter 10 of [4].Let us give the main ideas of the proof of Proposition 3.4.
(1) The sign of  `1{2 is a direct consequence of (3.13) and the choice of the initial condition  ´1{2 "  p 0 q: it is a convex combination of vectors with same signs.Moreover assuming  0 ě 0, the nonnegativity of   can be obtained by induction using the following arguments.The equation (3.14)Since ℎ{ 2 ă 1 and   ě 0, one has ˆ1 `ℎ 2  ˙ ě 0. Since   ě 0 and  `1{2 ě 0, one has that diag ``1{2 ˘ ě 0, so that the right-hand side of (3.17) is nonnegative componentwise.Moreover, the operator in the left-hand side of (3.17) has nonnegative inverse since it is an  -matrix for ℎ 1 P p0, ℎ 0 q small enough (depending on the bound on the maximum norm of the sequence p `1{2 q ě0 ).Indeed, one can check that it is a -matrix since its off-diagonal coefficients are nonpositive, and it is also a  -matrix (for ℎ P p0, ℎ 1 q).(2) Taking the scalar product of (3.14)  which implies the result (3.16).
We display in Figure 8 the comparison of the method above with the Lie splitting method, with the linear part approximated by a formula similar to (3.4), and with the implicit Euler method.We compute the  2 numerical errors using a reference solution obtained by a standard method of order 10 with a very small time step.Unsurprisingly, the three methods are of order 1 numerically.Moreover, the linearly implicit method is faster than the implicit Euler method for a given error, but slower than the Lie splitting method.Note that all the methods preserve the nonnegativity of the solution (as long as one has a bound on }  } 2 and ℎ is sufficiently small with respect to this bound, and under an additional CFL condition for the Lie splitting method).
We display in Figure 9 the plots of the initial datum and the final time solution obtained at  " 1 with the same linearly implicit method of order 1 (for ℎ " 1{p5 ˆ211 q) (left hand side) and the plot of the evolution of   (right hand side).This illustrates the results of Proposition 3.4: the numerical solution starting from a nonnegative initial datum stays nonnegative, and the discrete energy does not increase with time (see (3.16)).

Conclusion and perspectives
This paper introduces a new class of methods for the time integration of evolution problems set as systems of ODEs (or PDEs after space discretization).This class contains methods that are only linearly implicit, no matter the evolution equation.Moreover, the paper describes a specific way to design linearly implicit methods of any arbitrarily high order.Using suitable definitions of consistency and stability, we prove that such methods are actually of high order for ODEs, and the proof extends to finite systems of ODEs.We illustrate numerically that some of these methods are of the expected order for two examples of PDEs (nonlinear Schrödinger equation in 1d and 2d and a nonlinear heat equation in 1d), and discuss some of their qualitative properties.Our numerical results show that the linearly implicit methods introduced in this paper behave rather poorly in terms of efficiency for simple small systems of ODEs.In contrast, they illustrate numerically that a linearly implicit method of order 2 outperforms standard methods of order 2 from the litterature for a NLS equation on a domain where no spectral method can be applied.
Perspectives of this work include a rigorous analysis of these linearly implicit methods in PDE contexts (i.e.before discretization in space).In this direction, a recent result [7] proves that the relaxation method (2.5), which belongs to the class of methods presented here (see Rem. 2.1), is of order 2 when applied to the NLS equation.Another question is that of the possibility to design, in a systematic way, linearly implicit methods of high order with some qualitative properties adapted to the PDE problem (e.g.preservation of mass or energy for NLS equation, energy decrease for parabolic problems).

Figure
Figure Comparison of methods of order 1 applied to (3.1): On the left hand side: maximal numerical error as a function of the time step (logarithmic scales); on the right hand side: CPU time (in seconds) as a function of the maximal numerical error.

Figure
Figure Comparison of methods of order 2 applied to (3.1): On the left hand side: maximal numerical error as a function of the time step (logarithmic scales); on the right hand side: CPU time (in seconds) as a function of the maximal numerical error.

Figure 3 .
Figure 3.Comparison of methods of order 1 applied to (3.2): On the left hand side: maximal numerical error as a function of the time step (logarithmic scales); on the right hand side: CPU time (in seconds) as a function of the maximal numerical error.

Figure 4 .
Figure 4. Comparison of methods of order 2 applied to (3.2): On the left hand side: maximal numerical error as a function of the time step (logarithmic scales); on the right hand side: CPU time (in seconds) as a function of the maximal numerical error.

Figure 5 .
Figure 5.Comparison of methods of order 2 for mass and energy conservations: On the left hand side: variation of the mass with respect to time; on the right hand side: variation of the energy with respect to time.

Figure 6 .
Figure 6.2D domain used in the simulation of the nonlinear Schrödinger equation (3.7).

Figure 7 .
Figure 7.Comparison of methods of order 2 applied to (3.7): On the left hand side: maximal numerical error as a function of the time step (logarithmic scales); on the right hand side: CPU time (in seconds) as a function of the maximal numerical error.
The relation p ´ q  "     is equivalent to p1 ´ q  " p ` q  .Since    "  pq 1  1 where  pq 1is the first component of   , we infer that the relation p ´ q  "     is also equivalent to p1 ´ q  "  " 0, then, because the matrix  is strictly upper triangular and   ‰ 1, we have   " 0 and hence   is not an eigenvector.Therefore, if   is an eigenvector, Proof.Let us start with the estimate on  1  .First of all we use a Taylor expansion and write for all 1 ď  ď   pp  ` ℎqq "Let us denote by p  q the vector of R  with p ˝q p´1q p  q{p ´1q! as component number .
with  `1 ´ we obtain ´ } 2 2 " ´ p `1 ,  `1{2 q ` p  ,  ´1{2 q Figure 8.Comparison of methods of order 1 applied to (3.10): On the left hand side: maximal numerical error as a function of the time step (logarithmic scales); on the right hand side: CPU time (in seconds) as a function of the maximal numerical error.Figure 9. Initial datum and solution at the final time  " 1 with respect to space (left hand side) and evolution of   with respect to time (right hand side)." ´ ``1 ,  `1{2 ˘`  p  ,  ´1{2 q, (3.18)