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 timeintegration 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 semidiscrete 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 timeintegration methods preserve the above stability properties. Unless the timestep is sufficiently small, explicit timeintegration methods do not converge and/or give unphysical oscillations in the numerical solution. Even with sufficiently small timesteps, algebraic errors may increase unboundedly due to the large number of timesteps. The simplest example where the stability holds is the Euler implicit method,
where is the timestep. 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 firstorder accurate, i.e. its global time discretization error is . Therefore, to get a sufficiently small discretization error, one must choose very small timesteps, which means that the method becomes computationally expensive and also causes a stronger increase of roundoff errors. However, there exists stable timeintegration methods of arbitrary high order. They are of implicit RungeKutta quadrature type (see e.g.[15]), and belong to the class of Astable 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 RungeKutta 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,
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 twosided multiplication with , we get
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 timestep factor τ, but since we use higher order timediscretization methods, τ will not be very small and cannot damp out the inverse to some power of the spacediscretization parameter h that also occurs in .) Therefore, we choose . Note that this implies that .
The resulting relation (2) can now be written
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 .
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 secondorder diffusion or diffusionconvection 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 and V is a reflexive Banach space.
For proper functions and , then F is monotone, i.e.
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 timeintegration 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 Astable methods; see e.g.[3,7,8]. To get stability for all times and time steps, one requires where preferably . Such methods are called Lstable (Lambert) and stiffly Astable [3], respectively.
An important class of methods which are stiffly Astable is a particular class of the implicit RungaKutta 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
where is the solution at time t, , , , , and . The global discretization error of the component for this method is , i.e. it is a thirdorder method and it is stiffly Astable 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
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 secondorder factor can be preconditioned by the method in Section 2. The ability to factorize in secondorder factors and solve the arising systems as such twobytwo block matrix systems means that one only has to solve firstorder 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
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 [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
The spectral condition number is then bounded by
Proof Let be the product of . We have
It follows that
By the arithmeticgeometric means inequality, we have
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 RungeKutta 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 RungeKutta 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
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 EulerLagrange equations become then , that is,
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 EulerLagrange equation (12) is identical to the classical Newton’s law
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,
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
and
which equals the total energy of the system.
The corresponding EulerLagrange equations become now
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
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
where C is a symmetric matrix. For the solution of the Hamiltonian system (15), we shall use an implicit RungeKutta method based on Gauss quadrature.
The sstage RungeKutta method applied to an initial value problem, , is defined by
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
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 sstage 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
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,
If , then a Taylor expansion shows that
where takes values in a tube with radius about the solution .
It follows that if
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 timeintegration method with an adaptive time step control. We present here such a method based on the use of backward integration at each timestep using the adjoint operator. The use of adjoint operators in error estimates gives back to the classical AubinNitsche 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 Bstability, that holds for the numerical solution of implicit RungeKutta 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
Here, are the quadrature coefficients.
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
For the implicit RungeKutta method with approximate solution , it holds
where is a piecewise polynomial of degree m.
The corresponding residual equals
By the property of implicit RungeKutta 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
Hence,
Here, we can use the Galerkin orthogonality property (24) to get
where is a polynomial of degree m.
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 RungeKutta method, based on Gaussian quadrature, applied for hyperbolic (conservative) problems has order 2m.
6 A numerical test example
We consider the linear parabolic problem,
in the unit square domain with boundary condition
As initial function , we choose a tentlike 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 ninepoint 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 RungeKutta 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 timediscretized problems, we use the GMRES method with preconditioners from Section 2 and with the tolerance . Let us note that GMRES needs 56 iteration for this tolerance. The problem is implemented in Matlab.
The primary aim is to show how the timediscretization 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
In Figures 2 and 3, there are depicted numerical stationary solutions for and , respectively. The discretization parameter is .
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 timediscretization error decrease in dependence on τ at a fixed, relatively small, time . We consider five different timediscretization 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 timediscretization parameter , . So, we investigate the following error:
which values are found in Tables 35. 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 timediscretization 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τ
The error estimates from Tables 35 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 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 RungeKutta method with the twopoint 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 timediscretization error.
Table 7. Values of stabilized timeTin dependence onℓandτ
Table 8. Time discretization error at timein dependence onℓandτ
In Figures 4 and 5, there are depicted numerical stationary solutions for and , respectively. The discretization parameter is .
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 RungeKutta 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 RungeKutta methods. An alternative, stabilized explicit RungeKutta 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. (3531432/HiCi). The authors, therefore, acknowledge technical and financial support of KAU.
References

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

Butcher, JC: Implicit RungeKutta processes. Math. Comput.. 18, 50–64 (1964). Publisher Full Text

Axelsson, O: A class of Astable 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 Astable 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

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

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

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

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

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

Hairer, E, Lubich, Ch, Roche, M: The Numerical Solution of DifferentialAlgebraic Systems by RungeKutta Methods, Springer, Berlin (1989)

Petzold, LR: Order results for implicit RungeKutta 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 Astability. BIT. 16, 226–227 (1976). Publisher Full Text

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

Hundsdorfer, W, Verwer, JG: Numerical Solution of Time Dependent AdvectionDiffusionReaction Equations, Springer, Berlin (2003)