This article is part of the series Recent Advances in Operator Equations, Boundary Value Problems, Fixed Point Theory and Applications, and General Inequalities.

Open Access Research

On the solution of high order stable time integration methods

Owe Axelsson12, Radim Blaheta2, Stanislav Sysala2 and Bashir Ahmad1*

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


The electronic version of this article is the complete one and can be found online at: http://www.boundaryvalueproblems.com/content/2013/1/108


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

© 2013 Axelsson et al.; licensee Springer

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,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M1">View MathML</a>

arises. Here, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M2">View MathML</a>, M is a mass matrix and M, A are <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M3">View MathML</a> matrices. For a finite difference approximation, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M4">View MathML</a>, 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M5">View MathML</a> and converges to a fixed stationary solution as <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M6">View MathML</a>, independent of the initial value <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M7">View MathML</a>. 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,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M8">View MathML</a>

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M9">View MathML</a> is the time-step. Here, the eigenvalues of the inverse of the resulting matrix in the corresponding system,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M10">View MathML</a>

equal <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M11">View MathML</a> and satisfy the stability condition,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M12">View MathML</a>

Here, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M13">View MathML</a> 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M14">View MathML</a> for eigenvalues <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M15">View MathML</a>. This holds for the implicit Euler method, where

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M16">View MathML</a>

This method is only first-order accurate, i.e. its global time discretization error is <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M17">View MathML</a>. 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M18">View MathML</a> of the corresponding matrix B where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M19">View MathML</a>, and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M20">View MathML</a> is a linear function of <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M21">View MathML</a> at the quadrature points in the interval <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M22">View MathML</a>, satisfy <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M23">View MathML</a> for all normal matrices <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M24">View MathML</a> with <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M25">View MathML</a>. The highest order achieved, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M26">View MathML</a> occurs for Gauss quadrature where m equals to the number of quadrature points within each time interval.

To satisfy the second, desirable condition,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M27">View MathML</a>

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, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M28">View MathML</a>. 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M29">View MathML</a>. It shows also how the general case, where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M30">View MathML</a>, 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,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M31">View MathML</a>

(1)

We assume that M is spd and that <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M32">View MathML</a>, 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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M33">View MathML</a>

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M34">View MathML</a> is a parameter. We assume that A is a normal matrix, that is, has a full eigenvector space and further that the symmetric part, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M35">View MathML</a> of A is spd. To estimate the eigenvalues of <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M36">View MathML</a>, we write

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M37">View MathML</a>

After a two-sided multiplication with <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M38">View MathML</a>, we get

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M39">View MathML</a>

(2)

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M40">View MathML</a>, etc. and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M41">View MathML</a>. Note that, by similarity, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M42">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M43">View MathML</a> have the same eigenvalues.

We are interested in cases where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M44">View MathML</a> may have large eigenvalues. (In our application, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M44">View MathML</a> 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M46">View MathML</a>.) Therefore, we choose <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M47">View MathML</a>. Note that this implies that <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M48">View MathML</a>.

The resulting relation (2) can now be written

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M49">View MathML</a>

(3)

where

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M50">View MathML</a>

Since <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M51">View MathML</a>, the real part of the eigenvalues of <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M52">View MathML</a> are bounded above by 1. To find estimates of the eigenvalues <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M53">View MathML</a> of <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M54">View MathML</a>, let <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M55">View MathML</a> be eigensolutions of <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M44">View MathML</a>, i.e. let

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M57">View MathML</a>

It follows from (3) that for <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M58">View MathML</a>,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M59">View MathML</a>

We write <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M60">View MathML</a> so <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M61">View MathML</a>, where i is the imaginary unit. Note that <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M62">View MathML</a> so <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M63">View MathML</a>. Since, by assumption, the real part of μ is positive, it holds <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M64">View MathML</a>. A computation shows that the values of the factor <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M65">View MathML</a> are located in a disc in the complex plane with center at <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M66">View MathML</a> and radius <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M67">View MathML</a>, where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M68">View MathML</a>.

Hence, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M53">View MathML</a> is located in a disc with center at <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M70">View MathML</a> and radius <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M71">View MathML</a>.

For <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M72">View MathML</a>, i.e. for real eigenvalues of <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M44">View MathML</a>, then <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M74">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M75">View MathML</a>.

3 A stiffly stable time integration method

Consider a system of ordinary differential equations,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M76">View MathML</a>

(4)

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M77','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M77">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M78','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M78">View MathML</a>, M, A are <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M3">View MathML</a> 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M80','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M80">View MathML</a>, as <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M81','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M81">View MathML</a>.

Such stability results hold for more general problems, such as for a nonlinear parabolic problem,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M82','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M82">View MathML</a>

(5)

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M83','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M83">View MathML</a> and V is a reflexive Banach space.

For proper functions <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M84','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M84">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M85','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M85">View MathML</a>, then F is monotone, i.e.

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M86','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M86">View MathML</a>

(6)

Here, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M87','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M87">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M88">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M89','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M89">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M90','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M90">View MathML</a> denote the scalar product, and the corresponding norm in <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M91">View MathML</a>, respectively. In this case, one can easily derive the bound

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M92','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M92">View MathML</a>

where u, v are solution of (5) corresponding to different initial values. Consequently making use of the Gronwall lemma, we obtain

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M93','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M93">View MathML</a>

Hence, (5) is stable in this case.

If F is strongly monotone (or dissipative), i.e. (6) is valid with <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M94','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M94">View MathML</a>, then

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M95','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M95">View MathML</a>

i.e. (5) is asymptotically stable. In particular, the above holds for the test problem considered in Section 6.

For large eigenvalues of <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M96">View MathML</a>, 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M97','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M97">View MathML</a>, in (4), this amounts to proper approximations of the matrix exponential function <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M98','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M98">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M99','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M99">View MathML</a>, by a rational function,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M100','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M100">View MathML</a>

where

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M101','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M101">View MathML</a>

and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M102','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M102">View MathML</a> denotes eigenvalues by E. Furthermore, to cope with problems where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M103','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M103">View MathML</a>, but arbitrarily close to <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M104','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M104">View MathML</a>, one needs A-stable methods; see e.g.[3,7,8]. To get stability for all times and time steps, one requires <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M105','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M105">View MathML</a> where preferably <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M106','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M106">View MathML</a>. 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M107','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M107">View MathML</a>, where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M108','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M108">View MathML</a> are the Legendre polynomials, orthogonal on the interval <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M109','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M109">View MathML</a>, see e.g.[1] and references therein. Note that <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M110','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M110">View MathML</a> is a root for all <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M111','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M111">View MathML</a>. The case <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M112','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M112">View MathML</a> is identical to the implicit Euler method.

Following [5], we consider here the next simplest case, where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M113','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M113">View MathML</a>, for the numerical solution of (4) over a time interval <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M114','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M114">View MathML</a>.

In this case, the quadrature points (for a unit interval) are <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M115','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M115">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M116','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M116">View MathML</a> and the numerical solution <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M117','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M117">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M118','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M118">View MathML</a>, at <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M119','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M119">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M120','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M120">View MathML</a> satisfies

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M121','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M121">View MathML</a>

(7)

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M122','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M122">View MathML</a> is the solution at time t, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M123','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M123">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M124','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M124">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M125','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M125">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M126','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M126">View MathML</a>, and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M127','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M127">View MathML</a>. The global discretization error of the <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M128','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M128">View MathML</a>-component for this method is <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M129','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M129">View MathML</a>, i.e. it is a third-order method and it is stiffly A-stable even for arbitrary strong variations of the coefficient <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M130','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M130">View MathML</a>. 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M131','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M131">View MathML</a>, we derive a modified form of the Schur complement system, that involves only an inner system with matrix <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M132','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M132">View MathML</a>. To this end, but only for the derivation of the method, we scale first the system with the block diagonal matrix <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M133','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M133">View MathML</a> to get

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M134','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M134">View MathML</a>

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M135','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M135">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M136','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M136">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M137','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M137">View MathML</a>. The Schur complement system for <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M118','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M118">View MathML</a> is multiplied with <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M139','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M139">View MathML</a>. Using commutativity, we get then

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M140','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M140">View MathML</a>

or

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M141','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M141">View MathML</a>

Hence,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M142','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M142">View MathML</a>

where

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M143','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M143">View MathML</a>

(8)

For higher order Radau quadrature methods, the corresponding matrix polynomial in <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M144','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M144">View MathML</a> 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M145','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M145">View MathML</a> 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M146','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M146">View MathML</a> in the interval

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M147','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M147">View MathML</a>

(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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M148','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M148">View MathML</a>

(10)

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M34">View MathML</a> 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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M150','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M150">View MathML</a>

it holds

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M151','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M151">View MathML</a>

where

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M152','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M152">View MathML</a>

If<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M153','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M153">View MathML</a>, then<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M154','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M154">View MathML</a>and<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M155','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M155">View MathML</a>.

The spectral condition number is then bounded by

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M156','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M156">View MathML</a>

If<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M157','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M157">View MathML</a>, then

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M158','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M158">View MathML</a>

Proof Let <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M159','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M159">View MathML</a> be the <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M160','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M160">View MathML</a> product of <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M161','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M161">View MathML</a>. We have

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M162','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M162">View MathML</a>

It follows that

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M163','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M163">View MathML</a>

By the arithmetic-geometric means inequality, we have

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M164','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M164">View MathML</a>

(11)

a computation shows that

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M165','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M165">View MathML</a>

for <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M166','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M166">View MathML</a>, where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M167','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M167">View MathML</a>. Further, a computation shows that <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M168','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M168">View MathML</a>, which is in accordance with the lower bound in (11). Since

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M169','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M169">View MathML</a>

it follows that

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M170','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M170">View MathML</a>

or

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M171','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M171">View MathML</a>

For <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M172','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M172">View MathML</a>, a computation shows that

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M173','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M173">View MathML</a>

 □

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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M113','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M113">View MathML</a> 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M175','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M175">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M176','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M176">View MathML</a>. Then Radau quadrature on the interval <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M177','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M177">View MathML</a> has quadrature points <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M178','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M178">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M179','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M179">View MathML</a>, and coefficients <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M180','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M180">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M181','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M181">View MathML</a>, which results in the relation

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M182','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M182">View MathML</a>

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M183','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M183">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M184','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M184">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M185','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M185">View MathML</a> denote the corresponding approximations of u at <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M186','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M186">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M187','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M187">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M188','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M188">View MathML</a>, respectively.

This equation is coupled with an equation based on quadrature

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M189','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M189">View MathML</a>

which, using the stated quadrature rules, results in

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M190','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M190">View MathML</a>

that is,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M191','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M191">View MathML</a>

Remark 3.2 The arising system in a high order method involving <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M192','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M192">View MathML</a> 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, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M193','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M193">View MathML</a> as <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M194','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M194">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M195','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M195">View MathML</a>, where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M196','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M196">View MathML</a>, 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M197','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M197">View MathML</a>, 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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M198','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M198">View MathML</a>

with <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M199','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M199">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M200','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M200">View MathML</a> and it is normally assumed that the initial values satisfy the constraint equation, i.e.

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M201','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M201">View MathML</a>

If <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M202','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M202">View MathML</a> 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.

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M203','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M203">View MathML</a>

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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M204','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M204">View MathML</a>

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M205','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M205">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M206','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M206">View MathML</a>.

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,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M207','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M207">View MathML</a>

where K is the kinetic energy and V the potential energy. Here, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M208','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M208">View MathML</a> denote the Cartesian coordinate of the ith point mass <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M209','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M209">View MathML</a>.

The configuration strives to minimize the total energy. The corresponding Euler-Lagrange equations become then <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M210','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M210">View MathML</a>, that is,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M211','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M211">View MathML</a>

(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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M212','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M212">View MathML</a> according to

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M213','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M213">View MathML</a>

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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M214','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M214">View MathML</a>

Let <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M215','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M215">View MathML</a> be the momentum. Then

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M216','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M216">View MathML</a>

A mechanical system can be described by general coordinates

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M217','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M217">View MathML</a>

i.e. not necessarily Cartesian, but angles, length along a curve, etc. The Lagrangian takes the form <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M218','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M218">View MathML</a>. If q is determined to satisfy

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M219','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M219">View MathML</a>

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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M220','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M220">View MathML</a>

(13)

Letting here

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M221','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M221">View MathML</a>

be the momentum variable, and using the transformation <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M222','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M222">View MathML</a> we can write (13) as the Hamiltonian,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M223','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M223">View MathML</a>

For a mechanical system with potential energy a function of configuration only and kinetic energy K given by a quadratic form

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M224','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M224">View MathML</a>

where G is an spd matrix, possibly depending on q, we get

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M225','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M225">View MathML</a>

(14)

and

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M226','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M226">View MathML</a>

which equals the total energy of the system.

The corresponding Euler-Lagrange equations become now

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M227','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M227">View MathML</a>

(15)

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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M228','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M228">View MathML</a>

which, since <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M229','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M229">View MathML</a> implies <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M230','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M230">View MathML</a>, are hence equivalent to the Lagrange equations.

By (15), it holds

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M231','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M231">View MathML</a>

(16)

that is, the Hamiltonian function <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M232','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M232">View MathML</a> is a first integral for the system (15).

The flow <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M233','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M233">View MathML</a> of a Hamiltonian system is the mapping that describes the evolution of the solution by time, i.e.<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M234','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M234">View MathML</a>, where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M235','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M235">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M236','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M236">View MathML</a> is the solution of the system for the initial values <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M237','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M237">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M238','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M238">View MathML</a>.

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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M239','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M239">View MathML</a>

(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, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M240','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M240">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M241','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M241">View MathML</a> is defined by

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M242','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M242">View MathML</a>

(18)

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M243','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M243">View MathML</a>, see e.g.[1,4]. The familiar implicit midpoint rule is the special case where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M244','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M244">View MathML</a>. Here, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M245','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M245">View MathML</a> are the zeros of the shifted Legendre polynomial <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M246','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M246">View MathML</a>. 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M247','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M247">View MathML</a> is a polynomial of degree s, then (18) takes the form

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M248','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M248">View MathML</a>

(19)

and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M249','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M249">View MathML</a>.

For the Hamiltonian (17), it holds

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M250','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M250">View MathML</a>

and it follows from (16) that

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M251','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M251">View MathML</a>

Since the integrand is a polynomial of degree <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M252','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M252">View MathML</a>, it is evaluated exactly by the s-stage Gaussian quadrature formula. Therefore, since

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M253','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M253">View MathML</a>

it follows that the energy quadrature forms <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M254','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M254">View MathML</a> 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,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M255','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M255">View MathML</a>

holds, implying

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M256','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M256">View MathML</a>

(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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M257','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M257">View MathML</a>. Consider the problem <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M258','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M258">View MathML</a> where u belongs to some function space V and the corresponding truncation error,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M259','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M259">View MathML</a>

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M260','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M260">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M261','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M261">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M262','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M262">View MathML</a>.

If <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M263','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M263">View MathML</a>, then a Taylor expansion shows that

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M264','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M264">View MathML</a>

(21)

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M265','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M265">View MathML</a> takes values in a tube with radius <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M266','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M266">View MathML</a> about the solution <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M267','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M267">View MathML</a>.

It follows that if

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M268','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M268">View MathML</a>

(22)

and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M269','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M269">View MathML</a>, then

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M270','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M270">View MathML</a>

Under the above conditions, the discretization error <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M271','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M271">View MathML</a>, where

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M272','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M272">View MathML</a>

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M273','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M273">View MathML</a>, is the approximate solution, satisfies

(i) if F is strongly monotone and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M274','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M274">View MathML</a>, then <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M275','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M275">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M276','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M276">View MathML</a>;

(ii) if F is monotone (or conservative) and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M277','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M277">View MathML</a>, then <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M278','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M278">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M279','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M279">View MathML</a>.

Here, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M280','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M280">View MathML</a> depends on <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M281','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M281">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M282','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M282">View MathML</a>, but is independent of the stiffness of the problem under the appropriate conditions stated above.

If the solution u is smooth so that <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M283','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M283">View MathML</a> has also only smooth components, then <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M284','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M284">View MathML</a> may be much smaller than <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M285','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M285">View MathML</a>, showing that the stiffness, i.e. factors <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M286','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M286">View MathML</a>, do not enter in the error estimate.

In many problems, we can expect that <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M287','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M287">View MathML</a> is of the same order as <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M282','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M282">View MathML</a>, i.e. the first and last forms in (21) have the same order. In particular, this holds for a linear problem <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M289','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M289">View MathML</a>, where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M290','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M290">View MathML</a>.

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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M291','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M291">View MathML</a>-lifting method used in boundary value problems to derive discretization error estimates in <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M291','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M291">View MathML</a> 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M293','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M293">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M294','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M294">View MathML</a> be two approximate solutions to <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M295','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M295">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M279','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M279">View MathML</a> extended to polynomials of degree m from their pointwise values at <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M297','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M297">View MathML</a> in the interval <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M298','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M298">View MathML</a>. Let

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M299','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M299">View MathML</a>

Then, since by (18), <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M300','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M300">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M301','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M301">View MathML</a> satisfy the differential equation at the quadrature points, and by (19) it holds

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M302','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M302">View MathML</a>

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M303','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M303">View MathML</a>, where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M304','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M304">View MathML</a> is the set of quadrature points. Since <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M305','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M305">View MathML</a> is a polynomial of degree <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M306','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M306">View MathML</a>, Gauss quadrature is exact so

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M307','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M307">View MathML</a>

Here, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M308','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M308">View MathML</a> are the quadrature coefficients.

Hence,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M309','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M309">View MathML</a>

Since <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M310','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M310">View MathML</a>, 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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M311','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M311">View MathML</a>

(23)

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M312','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M312">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M313','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M313">View MathML</a>.

For the implicit Runge-Kutta method with approximate solution <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M314','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M314">View MathML</a>, it holds

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M315','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M315">View MathML</a>

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M316','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M316">View MathML</a> is a piecewise polynomial of degree m.

The corresponding residual equals

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M317','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M317">View MathML</a>

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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M318','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M318">View MathML</a>

(24)

to all polynomials of degree m. Here, the ‘dot’ indicates a vector product in <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M319','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M319">View MathML</a>. The discretization error equals <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M320','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M320">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M279','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M279">View MathML</a>. The error estimation will be based on the backward integration of the adjoint operator problem,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M322','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M322">View MathML</a>

(25)

Note that <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M323','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M323">View MathML</a>. It holds

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M324','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M324">View MathML</a>

so by integration by parts, we get

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M325','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M325">View MathML</a>

Here

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M326','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M326">View MathML</a>

Hence,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M327','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M327">View MathML</a>

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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M328','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M328">View MathML</a>

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M329','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M329">View MathML</a> is a polynomial of degree m.

Since <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M330','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M330">View MathML</a>, it follows that

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M331','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M331">View MathML</a>

and from <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M332','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M332">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M333','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M333">View MathML</a>, it follows that

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M334','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M334">View MathML</a>

Hence,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M335','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M335">View MathML</a>

Under sufficient regularity assumptions the last term can be bounded by <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M336','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M336">View MathML</a>. Hence, the discretization error grows linearly with time,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M337','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M337">View MathML</a>

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,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M338','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M338">View MathML</a>

(26)

in the unit square domain <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M339','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M339">View MathML</a> with boundary condition

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M340','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M340">View MathML</a>

(27)

As initial function <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M341','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M341">View MathML</a>, we choose a tent-like function with <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M342','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M342">View MathML</a> at the center of Ω and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M343','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M343">View MathML</a> on Ω; see Figure 1.

thumbnailFigure 1. Initial function.

Here, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M344','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M344">View MathML</a>, where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M345','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M345">View MathML</a> , <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M346','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M346">View MathML</a>, 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M347','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M347">View MathML</a> satisfies the conditions of the ratio <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M348','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M348">View MathML</a> from (9). We let <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M349','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M349">View MathML</a>.

Further b is a vector satisfying <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M350','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M350">View MathML</a>. We choose <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M351','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M351">View MathML</a>, where is a parameter, possibly <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M352','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M352">View MathML</a>.

After a finite element or finite difference approximation, a system of the form (4) arises. For a finite difference approximation <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M353','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M353">View MathML</a>, 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M354','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M354">View MathML</a> for <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M355','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M355">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M356','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M356">View MathML</a> for <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M357','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M357">View MathML</a>.

The time discretization is given by the implicit Runge-Kutta method with the Radau quadrature for <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M358','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M358">View MathML</a>; see Section 3. For comparison, we also consider <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M359','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M359">View MathML</a>, 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M360','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M360">View MathML</a>. 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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M361','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M361">View MathML</a>

(28)

6.1 Experiments with a known and smooth stationary solution

If we let

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M362','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M362">View MathML</a>

then the solution to (28) satisfies

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M363','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M363">View MathML</a>

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:

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M364','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M364">View MathML</a>

Here, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M365','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M365">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M366','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M366">View MathML</a> 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M367','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M367">View MathML</a>. 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M368','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M368">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M369','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M369">View MathML</a>, respectively. The discretization parameter is <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M370','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M370">View MathML</a>.

thumbnailFigure 2. Numerical stationary solution for<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M371','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M371">View MathML</a>.

thumbnailFigure 3. Numerical stationary solution for<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M372','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M372">View MathML</a>.

Now we will investigate how fast the numerical solution to (26)-(27) approaches the numerical solution to (28) in dependence on τ. We fix <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M373','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M373">View MathML</a> and we search the smallest time T for which

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M374','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M374">View MathML</a>

where the vectors <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M366','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M366">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M376','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M376">View MathML</a> 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 <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M378','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M378">View MathML</a>. We consider five different time-discretization parameters: <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M379','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M379">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M380','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M380">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M381','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M381">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M382','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M382">View MathML</a>, and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M383','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M383">View MathML</a>. We will compare the maximal differences between the vectors <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M384','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M384">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M385','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M385">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M386','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M386">View MathML</a>, where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M387','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M387">View MathML</a> represents the numerical solution to (26)-(27) at time T for the time-discretization parameter <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M388','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M388">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M389','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M389">View MathML</a>. So, we investigate the following error:

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M390','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M390">View MathML</a>

which values are found in Tables 3-5. If we let <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M373','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M373">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M369','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M369">View MathML</a> 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 time<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M378','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M378">View MathML</a>in dependence onhandτ

If we let <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M373','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M373">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M370','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M370">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M368','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M368">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M369','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M369">View MathML</a>, we obtain results in Table 4. We can see that the investigated time-discretization error decreases faster for <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M369','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M369">View MathML</a> than for <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M368','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M368">View MathML</a>.

Table 4. Time discretization error at time<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M378','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M378">View MathML</a>in dependence onandτ

If we let <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M401','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M401">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M373','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M373">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M370','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M370">View MathML</a>, and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M369','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M369">View MathML</a>, we obtain results in Table 5.

Table 5. Time discretization error at time<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M378','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M378">View MathML</a>in dependence onkandτ

The error estimates from Tables 3-5 indicate that the expected error estimate <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M406','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M406">View MathML</a> 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 time<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M378','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M378">View MathML</a>in dependence onandτfor the implicit Euler method

The error estimates are here significantly influenced by the oscillation parameter k. For the larger value <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M373','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M373">View MathML</a>, we do not observe convergence. In case <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M401','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M401">View MathML</a>, the convergence is first order <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M410','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M410">View MathML</a>, 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:

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M411','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M411">View MathML</a>

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 time<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M378','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M378">View MathML</a>in dependence onandτ

In Figures 4 and 5, there are depicted numerical stationary solutions for <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M368','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M368">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M414','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M414">View MathML</a>, respectively. The discretization parameter is <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M370','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M370">View MathML</a>.

thumbnailFigure 4. Numerical stationary solution for<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M371','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M371">View MathML</a>.

thumbnailFigure 5. Numerical stationary solution for<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M372','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/108/mathml/M372">View MathML</a>.

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 OpenURL

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

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

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

  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 OpenURL

  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 OpenURL

  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 OpenURL

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

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

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

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