Evolution equations arise in many important practical problems. They are frequently stiff, i.e. involves fast, mostly exponentially, decreasing and/or oscillating components. To handle such problems, one must use proper forms of implicit numerical time-integration methods. In this paper, we consider two methods of high order of accuracy, one for parabolic problems and the other for hyperbolic type of problems. For parabolic problems, it is shown how the solution rapidly approaches the stationary solution. It is also shown how the arising quadratic polynomial algebraic systems can be solved efficiently by iteration and use of a proper preconditioner.
Evolution equations arise in many important practical problems, such as for parabolic and hyperbolic partial differential equations. After application of a semi-discrete Galerkin finite element or a finite difference approximation method, a system of ordinary differential equations,
In the above applications, the order n of the system can be very large. Under reasonable assumptions of the given source function f, the system is stable, i.e. its solution is bounded for all and converges to a fixed stationary solution as , independent of the initial value . This holds if A is a normal matrix, that is, has a complete eigenvector space, and has eigenvalues with positive real parts. This condition holds for parabolic problems, where the eigenvalues of A are real and positive. In more involved problems, the matrix A may have complex eigenvalues with arbitrary large imaginary parts.
Clearly, not all numerical time-integration methods preserve the above stability properties. Unless the time-step is sufficiently small, explicit time-integration methods do not converge and/or give unphysical oscillations in the numerical solution. Even with sufficiently small time-steps, algebraic errors may increase unboundedly due to the large number of time-steps. The simplest example where the stability holds is the Euler implicit method,
Here, denotes the set of eigenvalues of A. To more quickly damp out initial transients in the solution, which arises for instance due to that the initial value may not satisfy boundary conditions given in the parabolic problem, one should preferably have eigenvalues of the inverse of the discrete matrix B, that satisfies for eigenvalues . This holds for the implicit Euler method, where
This method is only first-order accurate, i.e. its global time discretization error is . Therefore, to get a sufficiently small discretization error, one must choose very small time-steps, which means that the method becomes computationally expensive and also causes a stronger increase of round-off errors. However, there exists stable time-integration methods of arbitrary high order. They are of implicit Runge-Kutta quadrature type (see e.g.[1-5]), and belong to the class of A-stable methods, i.e. the eigenvalues of the corresponding matrix B where , and is a linear function of at the quadrature points in the interval , satisfy for all normal matrices with . The highest order achieved, occurs for Gauss quadrature where m equals to the number of quadrature points within each time interval.
To satisfy the second, desirable condition,
one can use a special subclass of such methods, based on Radau quadrature; see, e.g.[1,5]. The discretization error is here only one order less, . For linear problems, all such stable methods lead to rational polynomial approximation matrices B, and hence the need to solve quadratic polynomial equations. For stable methods, it turns out that the roots of these polynomials are complex.
In Section 2, a preconditioning method is described that is very efficient when solving such systems, without the need to factorize the quadratic polynomials in first order factors, thereby avoiding the need to use complex arithmetics. Section 3 discusses the special case where . It shows also how the general case, where , can be handled.
Section 4 deals with the use of implicit Runge-Kutta methods of Gauss quadrature type for solving hyperbolic systems of Hamiltonian type.
Section 5 presents a method to derive time discretization errors.
In Section 6, some illustrating numerical tests are shown. The paper ends with concluding remarks.
2 Preconditioners for quadratic matrix polynomials
From the introduction, it follows that it is of importance to use an efficient solution method for quadratic matrix polynomials and not factorize them in first order factors when this results in complex valued factors. For a method to solve complex valued systems in real arithmetics, see, e.g.. Here, we use a particular method that is suitable for the arising quadratic matrix polynomials.
Consider then the matrix polynomial,
We are interested in cases where may have large eigenvalues. (In our application, involves a time-step factor τ, but since we use higher order time-discretization methods, τ will not be very small and cannot damp out the inverse to some power of the space-discretization parameter h that also occurs in .) Therefore, we choose . Note that this implies that .
The resulting relation (2) can now be written
We write so , where i is the imaginary unit. Note that so . Since, by assumption, the real part of μ is positive, it holds . A computation shows that the values of the factor are located in a disc in the complex plane with center at and radius , where .
3 A stiffly stable time integration method
Consider a system of ordinary differential equations,
where , , M, A are matrices, where M is assumed to be spd and the symmetric part of A is positive semidefinite. In the practical applications that we consider, M corresponds to a mass matrix and A to a second-order diffusion or diffusion-convection matrix. Hence, n is large. Under reasonable assumptions on the source function f, such a system is stable for all t and its solution approaches a finite function, independent on the initial value , as .
Such stability results hold for more general problems, such as for a nonlinear parabolic problem,
where u, v are solution of (5) corresponding to different initial values. Consequently making use of the Gronwall lemma, we obtain
Hence, (5) is stable in this case.
i.e. (5) is asymptotically stable. In particular, the above holds for the test problem considered in Section 6.
For large eigenvalues of , such a system is stiff and can have fast decreasing and possibly oscillating components. This amounts to that the eigenvalues have large real part and possibly also large imaginary parts. To handle this, one needs stable numerical time-integration methods that do not contain corresponding increasing components. For , in (4), this amounts to proper approximations of the matrix exponential function , , by a rational function,
and denotes eigenvalues by E. Furthermore, to cope with problems where , but arbitrarily close to , one needs A-stable methods; see e.g.[3,7,8]. To get stability for all times and time steps, one requires where preferably . Such methods are called L-stable (Lambert) and stiffly A-stable , respectively.
An important class of methods which are stiffly A-stable is a particular class of the implicit Runga-Kutta methods; see [1,3,5]. Such methods correspond to rational polynomial approximations of the matrix exponential function with denominator having a higher degree than the nominator. Examples of such methods are based on Radau quadrature where the quadrature points are zeros of , where are the Legendre polynomials, orthogonal on the interval , see e.g. and references therein. Note that is a root for all . The case is identical to the implicit Euler method.
Following , we consider here the next simplest case, where , for the numerical solution of (4) over a time interval .
where is the solution at time t, , , , , and . The global discretization error of the -component for this method is , i.e. it is a third-order method and it is stiffly A-stable even for arbitrary strong variations of the coefficient . This can be compared with the trapezoidal or implicit midpoint methods which are only second order accurate and not stiffly stable.
The system in (7) can be solved via its Schur complement. Thereby, to avoid an inner system with matrix , we derive a modified form of the Schur complement system, that involves only an inner system with matrix . To this end, but only for the derivation of the method, we scale first the system with the block diagonal matrix to get
For higher order Radau quadrature methods, the corresponding matrix polynomial in is a mth order polynomial. By the fundamental theorem of algebra, one can factorize it in factors of at most second degree. They can be solved in a sequential order. Alternatively, using a method referred to in Remark 3.1, the solution components can be computed concurrently.
Each second-order factor can be preconditioned by the method in Section 2. The ability to factorize in second-order factors and solve the arising systems as such two-by-two block matrix systems means that one only has to solve first-order systems. This is of importance if for instance M and A are large sparse bandmatrices, since then one avoids increasing bandwidths in matrix products and one can solve systems of linear combinations of M and A more efficiently than for higher order polynomial combinations. Furthermore, this enables one to keep matrices on element by element form (see, e.g.) and it is in general not necessary to store the matrices M and A. The arising inner system can be solved by some inner iteration method.
Therefore, it is more efficient to keep the second order factors and instead solve the corresponding systems by preconditioned iterations. Thereby, the preconditioner involves only first order factors. As shown in Section 2, a very efficient preconditioner for the matrix B in (8) is
where is a parameter. As already shown in , for the above particular application it holds.
Proposition 3.1LetB, Cbe as defined in (8) and (10) and assume thatMis spd andAis spsd. Then letting
The spectral condition number is then bounded by
It follows that
By the arithmetic-geometric means inequality, we have
a computation shows that
it follows that
We conclude that the condition number is very close to its ideal unit value 1, leading to very few iterations. For instance, it suffices with at most 5 conjugate gradient iterations for a relative accuracy of 10−6.
For an early presentation of implicit Runge-Kutta methods, see  and also , where the method was called global integration method to indicate its capability for large values of m to use few, or even just one, time discretization steps. It was also shown that the coefficient matrix, formed by the quadrature coefficients had a dominating lower triangular part, enabling the use of a matrix splitting and Richardson iteration method. It can be of interest to point out that the Radau method for can be described in an alternative way, using Radau quadrature for the whole time step interval and combined with a trapezoidal method for the shorter interval.
This equation is coupled with an equation based on quadrature
which, using the stated quadrature rules, results in
Remark 3.2 The arising system in a high order method involving quadratic polynomial factors, can be solved sequentially in the order they appear. Alternatively (see, e.g., Exercise 2.31), one can use a method based on solving a matrix polynomial equation, as , , where , is the set of zeros of the polynomial and it is assumed that A has no eigenvalues in this set. (This holds in our applications.) Then, combining pairs of terms corresponding to complex conjugate roots , quadratic polynomials arise for the computation of the corresponding solution components. It is seen that in this method, the solution components can be computed concurrently.
Such a DAE is said to have index one, see e.g.. It can be seen to be a limit case of the system
Hence, such an DAE can be considered as an infinitely stiff differential equation problem. For strongly or infinitely stiff problems, there can occur an order reduction phenomenae. This follows since some high order error terms in the error expansion (cf. Section 5) are multiplied with (infinitely) large factors, leading to an order reduction for some methods. Heuristically, this can be understood to occur for the Gauss integration form of IRK but does not occur for the stiffly stable variants, such as based on the Radau quadrature. For further discussions of this, see, e.g.[10,13].
4 High order integration methods for Hamiltonian systems
Another important application of high order time integration methods occurs for Hamiltonian systems. Such systems occur in mechanics and particle physics, for instance. As an introduction, consider the conservation of energy principle. To this end, consider a mechanical system of k point masses and its associated Lagrangian functional,
This means that the Euler-Lagrange equation (12) is identical to the classical Newton’s law
A mechanical system can be described by general coordinates
then the motion of the system is described by the Lagrange equation,
For a mechanical system with potential energy a function of configuration only and kinetic energy K given by a quadratic form
where G is an spd matrix, possibly depending on q, we get
which equals the total energy of the system.
The corresponding Euler-Lagrange equations become now
and are referred to as the Hamiltonian system. This follows from
By (15), it holds
We consider now a Hamiltonian with a quadratic first integral in the form
where C is a symmetric matrix. For the solution of the Hamiltonian system (15), we shall use an implicit Runge-Kutta method based on Gauss quadrature.
where , see e.g.[1,4]. The familiar implicit midpoint rule is the special case where . Here, are the zeros of the shifted Legendre polynomial . For a linear problem, this results in a system which can be solved by the quadratic polynomial decomposition and the preconditioned iterative solution method, presented in Section 2.
For the Hamiltonian (17), it holds
and it follows from (16) that
This is an important property in Hamiltonian systems and is referred to as being symplectic. For further references of symplectic integrators, see .
5 Discretization error estimates
Error estimation methods for parabolic and hyperbolic problems can differ greatly. Parabolic problems are characterized by the monotonicity property (6) while for hyperbolic problems a corresponding conservation property,
Hence, there is no decrease of errors occurring at earlier time steps. On the other hand, the strong monotonicity property for parabolic problems implies that errors at earlier time steps decrease exponentially as time evolves.
For a derivation of discretization errors for such parabolic type problems for a convex combination of the implicit Euler method and the midpoint method, referred to as the θ-method, the following holds (see ). Similar estimates can also be derived for the Radau quadrature method, see, e.g..
The major result in  is the following.
It follows that if
It is seen from (20) that for hyperbolic (conservative) problems like the Hamiltonian problem in Section 4, the discretization error grows at least linearly with t, but likely faster if the solution is not sufficiently smooth. It may then be necessary to control the error by coupling the numerical time-integration method with an adaptive time step control. We present here such a method based on the use of backward integration at each time-step using the adjoint operator. The use of adjoint operators in error estimates gives back to the classical Aubin-Nitsche -lifting method used in boundary value problems to derive discretization error estimates in norm. It has also been used for error estimates in initial value problems, see e.g..
Assume that the monotonicity assumption (20) holds. We show first a nonlinear (monotone) stability property, called B-stability, that holds for the numerical solution of implicit Runge-Kutta methods based on Gauss quadrature points. It goes back to a scientific note in ; see also .
We present now a method for adaptive a posteriori error control for the initial value problem
The corresponding residual equals
By the property of implicit Runge-Kutta methods, it is orthogonal, i.e.
to all polynomials of degree m. Here, the ‘dot’ indicates a vector product in . The discretization error equals , . The error estimation will be based on the backward integration of the adjoint operator problem,
so by integration by parts, we get
Here, we can use the Galerkin orthogonality property (24) to get
i.e. the implicit Runge-Kutta method, based on Gaussian quadrature, applied for hyperbolic (conservative) problems has order 2m.
6 A numerical test example
We consider the linear parabolic problem,
As initial function , we choose a tent-like function with at the center of Ω and on ∂Ω; see Figure 1.
Figure 1. Initial function.
Here, , where , , is a parameter used to test the stability of the method with respect to oscillating coefficients. Here, τ is the time step to be used in the numerical solution of (26). Note that this function satisfies the conditions of the ratio from (9). We let .
After a finite element or finite difference approximation, a system of the form (4) arises. For a finite difference approximation , the identity matrix. The Laplacian operator is approximated with a nine-point difference scheme. We use an upwind discretization of the convection term. In the outer corner points of the domain, we use the boundary conditions for and for .
The time discretization is given by the implicit Runge-Kutta method with the Radau quadrature for ; see Section 3. For comparison, we also consider , i.e. the implicit Euler method, in some experiments. For solving the time-discretized problems, we use the GMRES method with preconditioners from Section 2 and with the tolerance . Let us note that GMRES needs 5-6 iteration for this tolerance. The problem is implemented in Matlab.
The primary aim is to show how the time-discretization errors decrease and how fast the numerical approximation of (26)-(27) approaches its stationary value, i.e. the corresponding numerical solution to the stationary problem
6.1 Experiments with a known and smooth stationary solution
If we let
then the solution to (28) satisfies
First, we will investigate the influence of the space discretization error on the stationary problem (28). To this end, we use the relative error estimate in the Euclidean norm:
Here, , denote the vectors representing the exact and numerical solutions to (28) at the nodal points, respectively. The error estimates in dependence on ℓ and h are found in Table 1. It is seen that the error decay is . This is caused by the use of first order upwind approximation of the convection term.
Table 1. The error estimates in dependence onℓandh
where the vectors and represent the numerical solution to (28) and the numerical solution to (26)-(27) at time T, respectively. The results for various ℓ and h are in Table 2. We can observe that the dependence of the results on h is small. For smaller ℓ, the final time does not depend on τ, while for larger ℓ, the dependence on τ is more significant.
Table 2. Values of timeTin dependence onhandτ
Finally, we investigate how the time-discretization error decrease in dependence on τ at a fixed, relatively small, time . We consider five different time-discretization parameters: , , , , and . We will compare the maximal differences between the vectors and , , where represents the numerical solution to (26)-(27) at time T for the time-discretization parameter , . So, we investigate the following error:
which values are found in Tables 3-5. If we let , and use various h, we obtain results written in Table 3. It is seen that the influence of h on the time discretization error is small for the larger time steps but more noticeable for the smaller time steps when the time and space discretization errors are of the same order.
Table 3. Time discretization error at timein dependence onhandτ
If we let , , and , we obtain results in Table 4. We can see that the investigated time-discretization error decreases faster for than for .
Table 4. Time discretization error at timein dependence onℓandτ
If we let , , , and , we obtain results in Table 5.
Table 5. Time discretization error at timein dependence onkandτ
Table 6. Time discretization error at timein dependence onℓandτfor the implicit Euler method
The error estimates are here significantly influenced by the oscillation parameter k. For the larger value , we do not observe convergence. In case , the convergence is first order , that is, much slower than for the Runge-Kutta method with the two-point Radau quadrature.
6.2 Experiments with an unknown and less smooth stationary solution
Here, we replace the above defined function g with the following one:
and prepare Tables 7 and 8 correspondingly to Tables 2 and 4, respectively. The results in Tables 7 and 8 are very similar to the results from Tables 2 and 4. It means that less smoothness in space of the solution to (26)-(27) do not significantly influence the time-discretization error.
Table 7. Values of stabilized timeTin dependence onℓandτ
Table 8. Time discretization error at timein dependence onℓandτ
7 Concluding remarks
There are several advantages in using high order time integration methods. Clearly, the major advantage is that the high order of discretization errors enables the use of larger, and hence fewer timesteps to achieve a desired level of accuracy. Some of the methods, like Radau integration, are highly stable, i.e. decrease unwanted solution components exponentially fast and do not suffer from an order reduction, which is otherwise common for many other methods. The disadvantage with such high order methods is that one must solve a number of quadratic matrix polynomial equations. For this reason, much work has been devoted to development of simpler methods, like diagonally implicit Runge-Kutta methods; see e.g.. Such methods are, however, of lower order and may suffer from order reduction.
In the present paper, it has been shown that the arising quadratic matrix system polynomial factors can be handled in parallel and each of them can be solved efficiently with a preconditioning method, resulting in very few iterations. Each iteration involves just two first order matrix real valued factors, similar to what arises in the diagonal implicit Runge-Kutta methods. An alternative, stabilized explicit Runge-Kutta methods, i.e. methods where the stability domain has been extended by use of certain forms of Chebyshev polynomials; see, e.g. can only be competitive for modestly stiff problems.
It has also been shown that the methods behave robustly with respect to oscillations in the coefficients in the differential operator. Hence, in practice, high order methods have a robust performance and do not suffer from any real disadvantage.
The authors declare that they have no competing interests.
Each of the authors, OA, RB, SS and BA, contributed to each part of this work equally and read and approved the final version of the manuscript.
This paper was funded by King Abdulaziz University, under grant No. (35-3-1432/HiCi). The authors, therefore, acknowledge technical and financial support of KAU.
Butcher, JC: Implicit Runge-Kutta processes. Math. Comput.. 18, 50–64 (1964). Publisher Full Text
Axelsson, O: A class of A-stable methods. BIT. 9, 185–199 (1969). Publisher Full Text
Axelsson, O: Global integration of differential equations through Lobatto quadrature. BIT. 4, 69–86 (1964). Publisher Full Text
Axelsson, O: On the efficiency of a class of A-stable methods. BIT. 14, 279–287 (1974). Publisher Full Text
Axelsson, O, Kucherov, A: Real valued iterative methods for solving complex symmetric linear systems. Numer. Linear Algebra Appl.. 7, 197–218 (2000). Publisher Full Text
Fried, I: Optimal gradient minimization scheme for finite element eigenproblems. J. Sound Vib.. 20, 333–342 (1972). Publisher Full Text
Petzold, LR: Order results for implicit Runge-Kutta methods applied to differential/algebraic systems. SIAM J. Numer. Anal.. 23(4), 837–852 (1986). Publisher Full Text
Axelsson, O: Error estimates over infinite intervals of some discretizations of evolution equations. BIT. 24, 413–424 (1984). Publisher Full Text
Wanner, G: A short proof on nonlinear A-stability. BIT. 16, 226–227 (1976). Publisher Full Text
Frank, R, Schneid, J, Ueberhuber, CW: The concept of B-convergence. SIAM J. Numer. Anal.. 18, 753–780 (1981). Publisher Full Text