# On the solution of high order stable time integration methods

Author Affiliations

1 King Abdulaziz University, Jeddah, Saudi Arabia

2 IT4 Innovations Department, Institute of Geonics AS CR, Ostrava, Czech Republic

For all author emails, please log on.

Boundary Value Problems 2013, 2013:108  doi:10.1186/1687-2770-2013-108

 Received: 15 February 2013 Accepted: 9 April 2013 Published: 26 April 2013

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### Abstract

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.

### 1 Introduction

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,

arises. Here, , M is a mass matrix and M, A are matrices. For a finite difference approximation, , the identity matrix.

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,

where is the time-step. Here, the eigenvalues of the inverse of the resulting matrix in the corresponding system,

equal and satisfy the stability condition,

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.[6]. Here, we use a particular method that is suitable for the arising quadratic matrix polynomials.

Consider then the matrix polynomial,

(1)

We assume that M is spd and that , which latter implies that the first order factors of B are complex. Systems with B will be solved by iteration. As a preconditioner, we use the matrix

where is a parameter. We assume that A is a normal matrix, that is, has a full eigenvector space and further that the symmetric part, of A is spd. To estimate the eigenvalues of , we write

After a two-sided multiplication with , we get

(2)

where , etc. and . Note that, by similarity, and have the same eigenvalues.

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

(3)

where

Since , the real part of the eigenvalues of are bounded above by 1. To find estimates of the eigenvalues of , let be eigensolutions of , i.e. let

It follows from (3) that for ,

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 .

Hence, is located in a disc with center at and radius .

For , i.e. for real eigenvalues of , then and .

### 3 A stiffly stable time integration method

Consider a system of ordinary differential equations,

(4)

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,

(5)

where and V is a reflexive Banach space.

For proper functions and , then F is monotone, i.e.

(6)

Here, , and , denote the scalar product, and the corresponding norm in , respectively. In this case, one can easily derive the bound

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.

If F is strongly monotone (or dissipative), i.e. (6) is valid with , then

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,

where

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 [3], 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.[1] and references therein. Note that is a root for all . The case is identical to the implicit Euler method.

Following [5], we consider here the next simplest case, where , for the numerical solution of (4) over a time interval .

In this case, the quadrature points (for a unit interval) are , and the numerical solution , , at and satisfies

(7)

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

where and , . The Schur complement system for is multiplied with . Using commutativity, we get then

or

Hence,

where

(8)

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.[9]) 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.

The problem with a direct factorization in first order factors is that complex matrix factors appear. This occurs for the matrix in (8) for a ratio of in the interval

(9)

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

(10)

where is a parameter. As already shown in [5], 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

it holds

where

If, thenand.

The spectral condition number is then bounded by

If, then

Proof Let be the product of . We have

It follows that

By the arithmetic-geometric means inequality, we have

(11)

a computation shows that

for , where . Further, a computation shows that , which is in accordance with the lower bound in (11). Since

it follows that

or

For , a computation shows 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.

Remark 3.1 High order implicit Runge-Kutta methods and their discretization error estimates can be derived using order tree methods as described in [1] and [10].

For an early presentation of implicit Runge-Kutta methods, see [2] and also [4], 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.

Namely, let , . Then Radau quadrature on the interval has quadrature points , , and coefficients , , which results in the relation

where , , denote the corresponding approximations of u at and and , respectively.

This equation is coupled with an equation based on quadrature

which, using the stated quadrature rules, results in

that is,

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.[11], 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.

Remark 3.3 Differential algebraic equations (DAE) arise in many important problems; see, for instance [10,12]. The simplest example of a DAE takes the form

with , and it is normally assumed that the initial values satisfy the constraint equation, i.e.

If in a sufficiently large set around the solution, one can formally eliminate the second part of the solution to form a differential equation in standard form.

Such a DAE is said to have index one, see e.g.[13]. It can be seen to be a limit case of the system

where and .

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,

where K is the kinetic energy and V the potential energy. Here, denote the Cartesian coordinate of the ith point mass .

The configuration strives to minimize the total energy. The corresponding Euler-Lagrange equations become then , that is,

(12)

We consider conservative systems, i.e. mechanical systems for which the total force on the elements of the system are related to the potential according to

This means that the Euler-Lagrange equation (12) is identical to the classical Newton’s law

Let be the momentum. Then

A mechanical system can be described by general coordinates

i.e. not necessarily Cartesian, but angles, length along a curve, etc. The Lagrangian takes the form . If q is determined to satisfy

then the motion of the system is described by the Lagrange equation,

(13)

Letting here

be the momentum variable, and using the transformation we can write (13) as the Hamiltonian,

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

(14)

and

which equals the total energy of the system.

The corresponding Euler-Lagrange equations become now

(15)

and are referred to as the Hamiltonian system. This follows from

which, since implies , are hence equivalent to the Lagrange equations.

By (15), it holds

(16)

that is, the Hamiltonian function is a first integral for the system (15).

The flow of a Hamiltonian system is the mapping that describes the evolution of the solution by time, i.e., where , is the solution of the system for the initial values , .

We consider now a Hamiltonian with a quadratic first integral in the form

(17)

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.

The s-stage Runge-Kutta method applied to an initial value problem, , is defined by

(18)

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.

If is a polynomial of degree s, then (18) takes the form

(19)

and .

For the Hamiltonian (17), it holds

and it follows from (16) that

Since the integrand is a polynomial of degree , it is evaluated exactly by the s-stage Gaussian quadrature formula. Therefore, since

it follows that the energy quadrature forms are conserved.

This is an important property in Hamiltonian systems and is referred to as being symplectic. For further references of symplectic integrators, see [10].

### 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,

holds, implying

(20)

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 [14]). Similar estimates can also be derived for the Radau quadrature method, see, e.g.[10].

The major result in [14] is the following.

Let . Consider the problem where u belongs to some function space V and the corresponding truncation error,

where , , .

If , then a Taylor expansion shows that

(21)

where takes values in a tube with radius about the solution .

It follows that if

(22)

and , then

Under the above conditions, the discretization error , where

, is the approximate solution, satisfies

(i) if F is strongly monotone and , then , ;

(ii) if F is monotone (or conservative) and , then , .

Here, depends on and , but is independent of the stiffness of the problem under the appropriate conditions stated above.

If the solution u is smooth so that has also only smooth components, then may be much smaller than , showing that the stiffness, i.e. factors , do not enter in the error estimate.

In many problems, we can expect that is of the same order as , i.e. the first and last forms in (21) have the same order. In particular, this holds for a linear problem , where .

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.[7].

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 [15]; see also [16].

Let , be two approximate solutions to , extended to polynomials of degree m from their pointwise values at in the interval . Let

Then, since by (18), and satisfy the differential equation at the quadrature points, and by (19) it holds

, where is the set of quadrature points. Since is a polynomial of degree , Gauss quadrature is exact so

Hence,

Since , this monotonicity property can be seen to hold also for the Radau quadrature method.

We present now a method for adaptive a posteriori error control for the initial value problem

(23)

where and .

For the implicit Runge-Kutta method with approximate solution , it holds

where is a piecewise polynomial of degree m.

The corresponding residual equals

By the property of implicit Runge-Kutta methods, it is orthogonal, i.e.

(24)

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,

(25)

Note that . It holds

so by integration by parts, we get

Here

Hence,

Here, we can use the Galerkin orthogonality property (24) to get

where is a polynomial of degree m.

Since , it follows that

and from and , it follows that

Hence,

Under sufficient regularity assumptions the last term can be bounded by . Hence, the discretization error grows linearly with time,

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,

(26)

in the unit square domain with boundary condition

(27)

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 .

Further b is a vector satisfying . We choose , where is a parameter, possibly .

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

(28)

#### 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 onandh

In Figures 2 and 3, there are depicted numerical stationary solutions for and , respectively. The discretization parameter is .

Figure 2. Numerical stationary solution for.

Figure 3. Numerical stationary solution for.

Now we will investigate how fast the numerical solution to (26)-(27) approaches the numerical solution to (28) in dependence on τ. We fix and we search the smallest time T for which

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 onandτ

If we let , , , and , we obtain results in Table 5.

Table 5. Time discretization error at timein dependence onkandτ

The error estimates from Tables 3-5 indicate that the expected error estimate holds.

For comparison, we perform the same experiment as in Table 5 for the implicit Euler time discretization. The results are in Table 6.

Table 6. Time discretization error at timein dependence onandτ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 onandτ

Table 8. Time discretization error at timein dependence onandτ

In Figures 4 and 5, there are depicted numerical stationary solutions for and , respectively. The discretization parameter is .

Figure 4. Numerical stationary solution for.

Figure 5. Numerical stationary solution for.

### 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.[10]. 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.[17] 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.

### Competing interests

The authors declare that they have no competing interests.

### Authors’ contributions

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.

### Acknowledgements

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.

### References

1. Butcher, JC: Numerical Method for Ordinary Differential Equations, Wiley, Chichester (2008)

2. Butcher, JC: Implicit Runge-Kutta processes. Math. Comput.. 18, 50–64 (1964). Publisher Full Text

3. Axelsson, O: A class of A-stable methods. BIT. 9, 185–199 (1969). Publisher Full Text

4. Axelsson, O: Global integration of differential equations through Lobatto quadrature. BIT. 4, 69–86 (1964). Publisher Full Text

5. Axelsson, O: On the efficiency of a class of A-stable methods. BIT. 14, 279–287 (1974). Publisher Full Text

6. 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

7. Varga, RS: Functional Analysis and Approximation Theory in Numerical Analysis, SIAM, Philadelphia (1971)

8. Gear, CW: Numerical Initial Value Problems in Ordinary Differential Equations, Prentice Hall, New York (1971)

9. Fried, I: Optimal gradient minimization scheme for finite element eigenproblems. J. Sound Vib.. 20, 333–342 (1972). Publisher Full Text

10. Hairer, E, Wanner, G: Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems, Springer, Berlin (1996)

11. Axelsson, O: Iterative Solution Methods, Cambridge University Press, Cambridge (1994)

12. Hairer, E, Lubich, Ch, Roche, M: The Numerical Solution of Differential-Algebraic Systems by Runge-Kutta Methods, Springer, Berlin (1989)

13. 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

14. Axelsson, O: Error estimates over infinite intervals of some discretizations of evolution equations. BIT. 24, 413–424 (1984). Publisher Full Text

15. Wanner, G: A short proof on nonlinear A-stability. BIT. 16, 226–227 (1976). Publisher Full Text

16. Frank, R, Schneid, J, Ueberhuber, CW: The concept of B-convergence. SIAM J. Numer. Anal.. 18, 753–780 (1981). Publisher Full Text

17. Hundsdorfer, W, Verwer, JG: Numerical Solution of Time Dependent Advection-Diffusion-Reaction Equations, Springer, Berlin (2003)