Multi-Adaptive Galerkin Methods for ODEs I

We present multi-adaptive versions of the standard continuous and discontinuous Galerkin methods for ODEs. Taking adaptivity one step further, we allow for individual time-steps, order and quadrature, so that in particular each individual component has its own time-step sequence. This paper contains a description of the methods, an analysis of their basic properties, and a posteriori error analysis. In the accompanying paper [A. Logg, SIAM J. Sci. Comput., 27 (2003), pp. 741-758], we present adaptive algorithms for time-stepping and global error control based on the results of the current paper.

1. Introduction. In this paper, we present multi-adaptive Galerkin methods for initial value problems for systems of ODEs of the form where u : [0, T ] → R N , f : R N × (0, T ] → R N is a given bounded function that is Lipschitz-continuous in u, u 0 ∈ R N is a given initial condition, and T > 0 is a given final time. We use the term multi-adaptivity to describe methods with individual time-stepping for the different components u i (t) of the solution vector u(t) = (u i (t)), including (i) time-step length, (ii) order, and (iii) quadrature, all chosen adaptively in a computational feedback process. In the companion paper [29], we apply the multi-adaptive methods to a variety of problems to illustrate the potential of multiadaptivity.
The ODE (1.1) models a very large class of problems, covering many areas of applications. Often different solution components have different time-scales and thus ask for individual time-steps. A prime example to be studied in detail below is our own solar system, where the moon orbits around Earth once every month, whereas the period of Pluto is 250 years. In numerical simulations of the solar system, the time-steps needed to track the orbit of the moon accurately are thus much less than those required for Pluto, the difference in time-scales being roughly a factor 3,000.
Surprisingly, individual time-stepping for ODEs has received little attention in the large literature on numerical methods for ODEs; see, e.g., [4,21,22,3,34]. For specific applications, such as the n-body problem, methods with individual time-stepping have been used-see, e.g., [31,1,5] or [25]-but a general methodology has been lacking. Our aim is to fill this gap. For time-dependent PDEs, in particular for conservation laws of the typeu + f (u) x = 0, attempts have been made to construct methods with individual (locally varying in space) time-steps. Flaherty et al. [20] have constructed a method based on the discontinuous Galerkin method combined with local forward Euler time-stepping. A similar approach is taken in [6], where a method based on the original work by Osher and Sanders [33] is presented for conservation laws in one and two space dimensions. Typically the time-steps used are based on local CFL conditions rather than error estimates for the global error and the methods are low order in time (meaning ≤ 2). We believe that our work on multi-adaptive Galerkin methods (including error estimation and arbitrary order methods) presents a general methodology to individual time-stepping, which will result in efficient integrators also for time-dependent PDEs.
The methods presented in this paper fall within the general framework of adaptive Galerkin methods based on piecewise polynomial approximation (finite element methods) for differential equations, including the continuous Galerkin method cG(q) of order 2q, and the discontinuous Galerkin method dG(q) of order 2q + 1; more precisely, we extend the cG(q) and dG(q) methods to their multi-adaptive analogues mcG(q) and mdG(q). Earlier work on adaptive error control for the cG(q) and dG(q) methods include [7,16,24,18,17,19]. The techniques for error analysis used in these references, developed by Johnson and coworkers (see, e.g., [11,12,10,13,14,15], and [8] in particular) naturally carries over to the multi-adaptive methods.
The outline of the paper is as follows: In section 2 we summarize the key features of the multi-adaptive methods, and in section 3 we discuss the benefits of the new methods in comparison to standard ODE codes. We then motivate and present the formulation of the multi-adaptive methods mcG(q) and mdG(q) in section 4. Basic properties of these methods, such as order, energy conservation, and monotonicity, are discussed in section 5. In the major part of this paper, section 6, we derive a posteriori error estimates for the two methods based on duality arguments, including Galerkin errors, numerical errors, and quadrature errors. We also prove an a posteriori error estimate for stability factors computed from approximate dual solutions.
2. Key features. We summarize the key features of our work on the mcG(q) and mdG(q) methods as follows.
2.1. Individual time-steps and order. To discretize (1.1), we introduce for each component, i = 1, . . . , N , a partition of the time-interval (0, T ] into M i subintervals, I ij = (t i,j−1 , t ij ], j = 1, . . . , M i , and we seek an approximate solution U (t) = (U i (t)) such that U i (t) is a polynomial of degree q ij on every local interval I ij . Each individual component U i (t) thus has its own sequence of time-steps, {k ij } Mi j=1 . The entire collection of individual time-intervals {I ij } may be organized into a sequence of time-slabs, collecting the time-intervals between certain synchronised time-levels common to all components, as illustrated in Figure 2.1.

Global error control.
Our goal is to compute an approximation U (T ) of the exact solution u(T ) at final time T within a given tolerance TOL > 0, using a minimal amount of computational work. This goal includes an aspect of reliability (the error should be less than the tolerance) and an aspect of efficiency (minimal computational work). To measure the error we choose a norm, such as the Euclidean norm · on R N , or more generally some other quantity of interest (see [32]).
The mathematical basis of global error control in · for mcG(q) is an error representation of the form solution U (t), ϕ(t) is the solution of an associated linearized dual problem, and (·, ·) is the R N scalar product.
Using the Galerkin orthogonality, the error representation can be converted into an error bound of the form are stability factors for the different components, depending on the dual solution ϕ(t), and where k i (t) = k ij , q i (t) = q ij for t ∈ I ij . The error bound may take different forms depending on how T 0 (R, ϕ) dt is bounded in terms of R and ϕ. By solving the dual problem numerically, the individual stability factors S i (T ) may be determined approximately, and thus the right-hand side of (2.2) may be evaluated. The adaptive algorithm seeks to satisfy the stopping criterion with maximal time-steps k = (k i (t)).

Iterative methods.
Both mcG(q) and mdG(q) give rise to systems of nonlinear algebraic equations, coupling the values of U (t) over each time-slab. Solving these systems with full Newton may be quite heavy, and we have instead successfully used diagonal Newton methods of more explicit nature.
2.4. Implementation of higher-order methods. We have implemented mcG(q) and mdG(q) in C++ for arbitrary q, which in practice means 2q ≤ 50. The implementation, Tanganyika, is described in more detail in [29] and is publicly (GNU GPL) available for Linux/Unix [30].
2.5. Applications. We have applied mcG(q) and mdG(q) to a variety of problems to illustrate their potential; see [29]. (See also [27] and [26].) In these applications, including the Lorenz system, the solar system, and a number of time-dependent PDE problems, we demonstrate the use of individual time-steps, and for each system we solve the dual problem to collect extensive information about the problems stability features, which can be used for global error control.
3. Comparison with standard ODE codes. Standard ODE codes use timesteps which are variable in time but the same for all components, and the time-steps are adaptively chosen by keeping the "local error" below a given local error tolerance set by the user. The global error connects to the local error through an estimate, corresponding to (2.2), of the form where S is a stability factor. Standard codes do not compute S, which means that the connection between the global error and the local error is left to be determined by the clever user, typically by computing with a couple of different tolerances.
Comparing the adaptive error control of standard ODE codes with the error control presented in this paper and the accompanying paper [29], an essential difference is thus the technique to estimate the global error: either by clever trial-and-error or, as we prefer, by solving the dual problem and computing the stability factors. Both approaches carry extra costs and what is best may be debated; see, e.g., [32] for a comparison.
However, expanding the scope to multi-adaptivity with individual stability factors for the different components, trial-and-error becomes very difficult or impossible, and the methods for adaptive time-stepping and error control presented below based on solving the dual problem seem to bring clear advantages in efficiency and reliability.
For a presentation of the traditional approach to error estimation in ODE codes, we refer to [2], where the following rather pessimistic view is presented: Here we just note that a precise error bound is often unknown and not really needed. We take the opposite view: global error control is always needed and often possible to obtain at a reasonable cost. We hope that multi-adaptivity will bring new life to the discussion on efficient and reliable error control for ODEs.

Multi-adaptive Galerkin.
In this section we present the multi-adaptive Galerkin methods, mcG(q) and mdG(q), based on the discretization presented in section 2.1.
4.1. The mcG(q) method. The mcG(q) method for (1.1) reads as follows: and where P q (I) denotes the linear space of polynomials of degree ≤ q on I. The trial functions in V are thus continuous piecewise polynomials, locally of degree q ij , and the test functions in W are discontinuous piecewise polynomials that are locally of degree q ij − 1.
Noting that the test functions are discontinuous, we can rewrite the global problem (4.1) as a number of successive local problems for each component: For i = 1, . . . , N , j = 1, . . . , M i , find U i | Iij ∈ P qij (I ij ) with U i (t i,j−1 ) given, such that We notice the presence of the vector U (t) = (U 1 (t), . . . , U N (t)) in the local problem for U i (t) on I ij . If thus component U i1 (t) couples to component U i2 (t) through f , this means that in order to solve the local problem for component U i1 (t) we need to know the values of component U i2 (t) and vice versa. The solution is thus implicitly defined by (4.3). Notice also that if we define the residual R of the approximate i.e., the residual is orthogonal to the test space on every local interval. We refer to this as the Galerkin orthogonality for the mcG(q) method.
Making an ansatz for every component U i (t) on every local interval I ij in terms of a nodal basis for P qij (I ij ) (see the appendix), we can rewrite (4.3) as m } q m=1 ⊂ P q−1 (0, 1) are corresponding polynomial weight functions, and τ ij maps Here we assume that the solution is expressed in terms of a nodal basis with the end-points included, so that by the continuity requirement ξ ij0 = ξ i,j−1,qi,j−1 .
Finally, evaluating the integral in (4.5) using nodal quadrature, we obtain a fully discrete scheme in the form of an implicit Runge-Kutta method: For i = 1, . . . , N , j = 1, . . . , M i , find {ξ ijm } qij m=0 , with ξ ij0 given by the continuity requirement, such that for certain weights {w n } (see the appendix). 4.2. The mdG(q) method. The mdG(q) method in local form, corresponding to (4.3), reads as follows: For i = 1, . . . , N , j = 1, . . . , M i , find U i | Iij ∈ P qij (I ij ), such that , and the initial condition is specified for i = 1, . . . , N , by U i (0 − ) = u i (0). On a global level, the trial and test spaces are given by In the same way as for the continuous method, we define the residual R of the approximate solution U as R i (U, t) =U i (t) − f i (U (t), t), defined on the inner of every local interval I ij , and we rewrite (4.7) in the form We refer to this as the Galerkin orthogonality for the mdG(q) method. Notice that this is similar to (4.4) if we extend the integral in (4.4) to include the left end-point of the interval I ij . (The derivative of the discontinuous solution is a Dirac delta function at the end-point.) Making an ansatz for the solution in terms of some nodal basis, we get, as for the continuous method, the following explicit version of (4.7) on every local interval: or, applying nodal quadrature, where the weight functions, the nodal points, and the weights are not the same as for the continuous method.

4.3.
The multi-adaptive mcG(q)-mdG(q) method. The discussion above for the two methods extends naturally to using different methods for different components. Some of the components could therefore be solved for using the mcG(q) method, while for others we use the mdG(q) method. We can even change methods between different intervals.
Although the formulation thus includes adaptive orders and methods, as well as adaptive time-steps, our focus will be mainly on adaptive time-steps.

Choosing basis functions and quadrature.
What remains in order to implement the two methods specified by (4.6) and (4.11) is to choose basis functions and quadrature. For simplicity and efficiency reasons, it is desirable to let the nodal points for the nodal basis coincide with the quadrature points. It turns out that for both methods, the mcG(q) and the mdG(q) methods, this is possible to achieve in a natural way. We thus choose the q + 1 Lobatto quadrature points for the mcG(q) method, i.e., the zeros of xP q (x) − P q−1 (x), where P q is the qth-order Legendre polynomial on the interval; for the mdG(q) method, we choose the Radau quadrature points, i.e., the zeros of P q (x) + P q+1 (x) on the interval (with time reversed so as to include the right end-point). See [28] for a detailed discussion on this subject. The resulting discrete schemes are related to the implicit Runge-Kutta methods referred to as Lobatto and Radau methods; see, e.g., [3].
5. Basic properties of the multi-adaptive Galerkin methods. In this section we examine some basic properties of the multi-adaptive methods, including order, energy conservation, and monotonicity.

5.1.
Order. The standard cG(q) and dG(q) methods are of order 2q and 2q + 1, respectively. The corresponding properties hold for the multi-adaptive methods, i.e., mcG(q) is of order 2q and mdG(q) is of order 2q + 1, assuming that the exact solution u is smooth. We examine this in more detail in subsequent papers.

Energy conservation for mcG(q).
The standard cG(q) method is energy-conserving for Hamiltonian systems. We now prove that also the mcG(q) method has this property, with the natural restriction that we should use the same time-steps for every pair of positions and velocities. We consider a Hamiltonian system, together with initial conditions for x andẋ. Hereẍ is the acceleration, which by Newton's second law is balanced by the force F (x) = −∇ x P (x) for some potential field P . With u = x and v =ẋ we rewrite (5.1) as The total energy E(t) is the sum of the kinetic energy K(t) and the potential energy P (x(t)), Multiplying (5.1) withẋ it is easy to see that energy is conserved for the continuous problem, i.e., E(t) = E(0) for all t ∈ [0, T ]. We now prove the corresponding property for the discrete mcG(q) solution of (5.2). Theorem 5.1. The multi-adaptive continuous Galerkin method conserves energy in the following sense: Let (U, V ) be the mcG(q) solution to (5.2) defined by (4.3). Assume that the same time-steps are used for every pair of positions and corresponding velocities. Then at every synchronized time-levelt, such as, e.g., T , we have 2 and P (t) = P (U (t)). Proof. If every pair of positions and velocities have the same time-step sequence, then we may chooseV as a test function in the equations for U , to get Similarly,U may be chosen as a test function in the equations for V to get and thus K(t) + P (t) = K(0) + P (0). Remark 5.1. Energy conservation requires exact integration of the right-hand side f (or at least that t 0 (U ,V ) dt + (P (t) − P (0)) = 0) but can also be obtained in the case of quadrature; see [23].

Monotonicity.
We shall prove that the mdG(q) method is B-stable (see [3]).
Theorem 5.2. Let U and V be the mdG(q) solutions of (1.1) with initial data U (0 − ) and V (0 − ), respectively, defined by (4.7) on the same partition. If the righthand side f is monotone, i.e., then, at every synchronized time-levelt, such as, e.g., T , we have Proof. Choosing the test function as v = W = U − V in (4.7) for U and V , summing over the local intervals, and subtracting the two equations, we have Noting that The proof is completed noting that the same analysis applies with T replaced by any other synchronized time-levelt. Remark 5.2. The proof extends to the fully discrete scheme, using the positivity of the quadrature weights.

6.
A posteriori error analysis. In this section we prove a posteriori error estimates for the multi-adaptive Galerkin methods, including quadrature and discrete solution errors. Following the procedure outlined in the introduction, we first define the dual linearized problem and then derive a representation formula for the error in terms of the dual and the residual.
6.1. The dual problem. The dual problem comes in two different forms: a continuous and a discrete. For the a posteriori error analysis of this section, we will make use of the continuous dual. The discrete dual problem is used to prove a priori error estimates.
To set up the continuous dual problem, we define, for given functions v 1 (t) and v 2 (t), where * denotes the transpose, and we note that The continuous dual problem is then defined as the following system of ODEs: with data ϕ T and right-hand side g. Choosing the data and right-hand side appropriately, we obtain error estimates for different quantities of the computed solution. We shall assume below that the dual solution has q continuous derivatives (ϕ (qij ) i ∈ C(I ij ) locally on interval I ij ) for the continuous method and q + 1 continuous derivatives (ϕ (qij +1) i ∈ C(I ij ) locally on interval I ij ) for the discontinuous method.

Error representation.
The basis for the error analysis is the following error representation, expressing the error of an approximate solution U (t) in terms of the residual R(U, t) via the dual solution ϕ(t). We stress that the result of the theorem is valid for any piecewise polynomial approximation of the solution to the initial value problem (1.1) and thus in particular the mcG(q) and mdG(q) approximations.
Theorem 6.1. Let U be a piecewise polynomial approximation of the exact solution u of (1.1), and let ϕ be the solution to (6.3) with right-hand side g(t) and initial data ϕ T , and define the residual of the approximate solution U as R( Assume also that U is right-continuous at T . Then the error e = U − u satisfies (6.4) Proof. By the definition of the dual problem, we have using (6.2) Integrating by parts, we get Thus, with L ϕT ,g (e) = (e(T ), ϕ T ) + T 0 (e, g) dt, we have which completes the proof. We now apply this theorem to represent the error in various norms. As before, we let · denote the Euclidean norm on R N and define v L 1 ([0,T ],R n ) = T 0 v dt. Corollary 6.2. If ϕ T = e(T )/ e(T ) and g = 0, then Corollary 6.3. If ϕ T = 0 and g(t) = e(t)/ e(t) , then

Galerkin errors.
To obtain expressions for the Galerkin errors, i.e., the errors of the mcG(q) or mdG(q) approximations, assuming exact quadrature and exact solution of the discrete equations, we use two ingredients: the error representation of Theorem 6.1 and the Galerkin orthogonalities, (4.4) and (4.9). We first prove the following interpolation estimate.
, then there is a constant C q , depending only on q, such that where π [q] f (x) is the qth-order Taylor expansion of f around x 0 = (a+b)/2, k = b−a, and C q = 1/(2 q q!).
Proof. Using Taylor's formula with the remainder in integral form, we have Note that since we allow the polynomial degree to change between different components and between different intervals, the interpolation constant will change in the same way. We thus have C qi = C qi (t) = C qij for t ∈ I ij .
We can now prove a posteriori error estimates for the mcG(q) and mdG(q) methods. The estimates come in a number of different versions. We typically use E 2 or E 3 to adaptively determine the time-steps and E 0 or E 1 to evaluate the error. The quantities E 4 and E 5 may be used for qualitative estimates of error growth. We emphasize that all of the estimates derived in Theorems 6.5 and 6.6 below may be of use in an actual implementation, ranging from the very sharp estimate E 0 containing only local quantities to the more robust estimate E 5 containing only global quantities.
Theorem 6.5. The mcG(q) method satisfies the following estimates: where (6.10) with C q as in Lemma 6.4, k i (t) = k ij , r i (t) = r ij , and s for t ∈ I ij , (6.11) and where π k ϕ is any test space approximation of the dual solution ϕ. Expressions such as C q−1 k q R are defined componentwise, i.e., Proof. Using the error representation of Theorem 6.1 and the Galerkin orthogonality (4.4), noting that the jump terms disappear since U is continuous, we have where π k ϕ is any test space approximation of ϕ. By the triangle inequality, we have Choosing π k ϕ i as in Lemma 6.4 on every interval I ij , we have Continuing, we have , and, finally, As an alternative we can use Cauchy's inequality in a different way. Continuing from E 2 , we have where |R(U, ·)| denotes the vector-valued function with components |R| i = |R i | = |R i (U, ·)|. Noting now that s is the L 2 -projection of |ϕ (q) | onto the piecewise constants on the partition, we have The proof of the estimates for the mdG(q) method is obtained similarly. Since in the discontinuous method the test functions are on every interval of one degree higher order than the test functions in the continuous method, we can choose a better interpolant. Thus, in view of Lemma 6.4, we obtain an extra factor k ij in the error estimates.
Theorem 6.6. The mdG(q) method satisfies the following estimates: where (6.14) and we otherwise use the notation of Theorem 6.5. Proof. As in the proof for the continuous method, we use the error representation of Theorem 6.1 and the Galerkin orthogonality (4.9) to get By Lemma 6.4 we obtain Continuing now in the same way as for the continuous method, we have E 2 ≤ E 3 ≤ E 4 and E 2 ≤ E 5 . Remark 6.1. When evaluating the expressions E 0 or E 1 , the interpolant π k ϕ does not have to be chosen as in Lemma 6.4. This is only a convenient way to obtain the interpolation constant. In section 6.6 below we discuss a more convenient choice of interpolant.
Remark 6.2. If we replace 1 kij Iij |R i | dt by max Iij |R i |, we may replace C q by a smaller constant C ′ q . The value of the constant thus depends on the specific way the residual is measured.

Computational errors.
The error estimates of Theorems 6.5 and 6.6 are based on the Galerkin orthogonalities (4.4) and (4.9). If the corresponding discrete equations are not solved exactly, there will be an additional contribution to the total error. Although Theorem 6.1 is still valid, the first steps in Theorems 6.5 and 6.6 are not. Focusing on the continuous method, the first step in the proof of Theorem 6.5 is the subtraction of a test space interpolant. This is possible, since by the Galerkin orthogonality we have for all test space interpolants π k ϕ. If the residual is no longer orthogonal to the test space, we add and subtract this term to get to the point where the implications of Theorem 6.5 are valid for one of the terms. Assuming now that ϕ varies slowly on each subinterval, we estimate the remaining extra term as follows: (6.16) whereφ is a piecewise constant approximation of ϕ (using, say, the mean values on the local intervals), is a stability factor, and we define the discrete or computational residual as More precise estimates may be used if needed. For the mdG(q) method, the situation is similar with the computational residual now defined as Thus, to estimate the computational error, we evaluate the computational residuals and multiply with the computed stability factors. 6.5. Quadrature errors. We now extend our analysis to take into account also quadrature errors. We denote integrals evaluated by quadrature with˜ . Starting from the error representation as before, we have for the mcG(q) method (6.20) if the quadrature is exact forU v when v is a test function. The first term of this expression was estimated in Theorem 6.5 and the second term is the computational error discussed previously (where˜ denotes that in a real implementation, (6.18) is evaluated using quadrature). The third term is the quadrature error, which may be nonzero even if f is linear, if the time-steps are different for different components. To estimate the quadrature error, notice that (6.21) are the same stability factors as in the estimate for the computational error and is the quadrature residual. A similar estimate holds for the mdG(q) method. We now make a few comments on how to estimate the quadrature residual. The Lobatto quadrature of the mcG(q) method is exact for polynomials of degree less than or equal to 2q − 1, and we have an order 2q estimate for˜ − in terms of f (2q) , and so we make the assumption R Q ij ∝ k 2qij ij . If, instead of using the standard quadrature rule over the interval with quadrature residual R Q0 ij , we divide the interval into 2 m parts and use the quadrature on every interval, summing up the result, we will get a different quadrature residual, namely where we have dropped the ij subindices. Thus, since Thus, by computing the integrals at two or more dyadic levels, we may estimate quadrature residuals and thus the quadrature error.
For the mdG(q) method the only difference is that the basic quadrature rule is one order better, i.e., instead of 2q we have 2q + 1, so that 6.6. Evaluating E G . We now present an approach to estimating the quantity T 0 (R(U, ·), ϕ − π k ϕ) dt by direct evaluation, with ϕ a computed dual solution and π k ϕ a suitably chosen interpolant. In this way we avoid introducing interpolation constants and computing derivatives of the dual. Note, however, that although we do not explicitly compute any derivatives of the dual, the regularity assumed in section 6.1 for the dual solution is still implicitly required for the computed quantities to make sense. Starting now with for the continuous method, we realize that the best possible choice of interpolant, if we want to prevent cancellation, is to choose π k ϕ such that R i (U, ·)(ϕ i − π k ϕ i ) ≥ 0 (or ≤ 0) on every local interval I ij . With such a choice of interpolant, we would have (6.27) The following lemmas give us an idea of how to choose the interpolant. Lemma 6.7. If, for i = 1, . . . , N , , t) and f i is linear or, alternatively, f = f (U (t), t) is linear and all components have the same timesteps and order, then every component R i (U, ·) of the mcG(q) residual is a Legendre polynomial of order q ij on I ij , for j = 1, . . . , M i .
Proof. On every interval I ij the residual component R i (U, ·) is orthogonal to P qij −1 (I ij ). Since the conditions assumed in the statement of the lemma guarantee that the residual is a polynomial of degree q ij on every interval I ij , it is clear that on every such interval it is the q ij th-order Legendre polynomial (or a multiple thereof).
Even if the rather strict conditions of this lemma do not hold, we can say something similar. The following lemma restates this property in terms of approximations of the residual. Lemma 6.8. LetR be the local L 2 -projection of the mcG(q) residual R onto the trial space, i.e.,R i (U, ·)| Iij is the L 2 (I ij )-projection onto P qij (I ij ) of R i (U, ·)| Iij , j = 1, . . . , M i , i = 1, . . . , N . Then everyR i (U, ·)| Iij is a Legendre polynomial of degree q ij . Proof.
for all v ∈ P qij −1 (I ij ), so thatR i (U, ·) is the q ij th-order Legendre polynomial on I ij .
Proof. We can write the qth-order Legendre polynomial on [−1, 1] as P q (x) = 1 q!2 q D q ((x 2 − 1) q ). Thus, integrating by parts, we have is zero at −1 and 1 for l < q. Assume now that f is a polynomial of degree q on [−1, 1] with the property (6.28). Since are linearly independent on [−1, 1] and orthogonal to the Radau polynomial Q q , {Q q (x), (x + 1), (x + 1) 2 , . . . , (x + 1) q } form a basis for P q ([−1, 1]). If then f is orthogonal to the subspace spanned by {(x + 1) p } q p=1 , we must have f = cQ q for some constant c, and the proof is complete.
, t) and f i is linear or, alternatively, f = f (U (t), t) is linear and all components have the same timesteps and order, then every component R i (U, ·) of the mdG(q) residual is a Radau polynomial of order q ij on I ij for j = 1, . . . , M i .
Proof. Note first that by assumption the residual R i (U, ·) is a polynomial of degree q ij on I ij . By the Galerkin orthogonality, we have which holds especially for v(t) = (t − t i,j−1 ) p with p = 1, . . . , q, for which the jump terms disappear. Rescaling to [−1, 1], it follows from Lemma 6.9 that the residual R i (U, ·) must be a Radau polynomial on I ij . Also for the discontinuous method there is a reformulation in terms of approximations of the residual. Lemma 6.11. LetR be the local L 2 -projection of the mdG(q) residual R onto the trial space, i.e.,R i (U, ·)| Iij is the L 2 (I ij )-projection onto P qij (I ij ) of R i (U, ·)| Iij , j = 1, . . . , M i , i = 1, . . . , N . Then everyR i (U, ·)| Iij is a Radau polynomial of degree q ij .
Proof. SinceR i (U, ·) is the L 2 -projection of R i (U, ·) onto P qij (I ij ) on I ij , it follows from the Galerkin orthogonality that for any v(t) = (t − t i,j−1 ) p with 1 ≤ p ≤ q. From Lemma 6.9 it then follows that R i (U, ·) is a Radau polynomial on I ij .
We thus know that the mcG(q) residuals are (in the sense of Lemma 6.8) Legendre polynomials on the local intervals and that the mdG(q) residuals are (in the sense of Lemma 6.11) Radau polynomials. This is illustrated in Figure 6.1.
From this information about the residual, we now choose the interpolant. Assume that the polynomial order of the method on some interval is q for the continuous method. Then the dual should be interpolated by a polynomial of degree q −1, i.e., we have freedom to interpolate at exactly q points. Since a qth-order Legendre polynomial has q zeros on the interval, we may choose to interpolate the dual exactly at those points where the residual is zero. This means that if the dual can be approximated well enough by a polynomial of degree q, the product R i (U, ·)(ϕ i − π k ϕ i ) does not change sign on the interval.
For the discontinuous method, we should interpolate the dual with a polynomial of degree q, i.e., we have freedom to interpolate at exactly q + 1 points. To get rid of the jump terms that are present in the error representation for the discontinuous method, we want to interpolate the dual at the beginning of every interval. This leaves q degrees of freedom. We then choose to interpolate the dual at the q points within the interval where the Radau polynomial is zero.
As a result, we may choose the interpolant in such a way that we have (6.29) with α ij = ±1, for both the mcG(q) method and the mdG(q) method (but the interpolants are different). Notice that the jump terms for the discontinuous method have disappeared.
There is now a simple way to compute the integrals Iij R i (U, ·)(ϕ i − π k ϕ i ) dt. Since the integrands are, in principle, products of two polynomials for which we know the positions of the zeros, the product is a polynomial with known properties. There are then constants C q (which can be computed numerically), depending on the order and the method, such that (6.30) Finally, note that there are "computational" counterparts also for the estimates of type E 3 in Theorems 6.5 and 6.6, namely for the continuous method and similarly for the discontinuous method. 6.7. The total error. The total error is composed of three parts-the Galerkin error, E G , the computational error, E C and the quadrature error, E Q : As an example, choosing estimate E 3 of Theorems 6.5 and 6.6 we have the following (approximate) error estimate for the mcG(q) method: for the mdG(q) method we have (6.34) These estimates containing Galerkin errors, computational errors, and quadrature errors also include numerical round-off errors (included in the computational error).
Modelling errors could also be similarly accounted for since these are closely related to quadrature errors, in that both errors can be seen as arising from integrating the wrong right-hand side.
The true global error may thus be estimated in terms of computable stability factors and residuals. We expect the estimate for the Galerkin error, E G , to be quite sharp, while E C and E Q may be less sharp. Even sharper estimates are obtained using estimates E 0 , E 1 , or E 2 of Theorems 6.5 and 6.6. 6.8. An a posteriori error estimate for the dual. We conclude this section by proving a computable a posteriori error estimate for the dual problem. To compute the stability factors used in the error estimates presented above, we solve the dual problem numerically, and we thus face the problem of estimating the error in the stability factors.
To demonstrate how relative errors of stability factors can be estimated using the same technique as above, we compute the relative error for the stability factor S ϕ (T ), defined as for a computed approximation Φ of the dual solution ϕ.
To estimate the relative error of the stability factor, we use the error representation of Theorem 6.1 to represent the L 1 ([0, T ], R N )-error of Φ in terms of the residual of Φ and the dual of the dual, ω. In [28] we prove the following lemma, from which the estimate follows.
Lemma 6.12. Let ϕ be the dual solution with stability factor S ϕ (t), i.e., with data ϕ(t) = 1 specified at time t, and let ω be the dual of the dual. We then have the following estimate: Theorem 6.13. Let Φ be a continuous approximation of the dual solution with residual R Φ , and assume that S ϕ (t)/S ϕ (T ) is bounded by C on [0, T ]. Then the following estimate holds for the relative error of the stability factor S Φ (T ): and for many problems we may take C = 1.
Proof. By Corollary 6.3, we have an expression for the L 1 ([0, T ], R N )-error of the dual, so that (6.38) With C defined as above it now follows by Lemma 6.12 that and the proof is complete. Remark 6.3. We also have to take into account quadrature errors when evaluating (6.35). This can be done in many ways; see, e.g., [9].
Appendix A. Derivation of the methods. This section contains some details left out of the discussion of section 4.
A.1. The mcG(q) method. To rewrite the local problem in a more explicit form, let {s n } q n=0 be a set of nodal points on [0, 1], with s 0 = 0 and s q = 1. A good choice for the cG(q) method is the Lobatto points of [0, 1]. Now, let τ ij be the linear mapping from the interval I ij to (0, 1], defined by n } q n=0 be the {s n } q n=0 Lagrange basis functions for P q ([0, 1]) on [0, 1], i.e., .
We can then express U i on I ij in the form and choosing the λ To simplify the notation, we drop the ij subindices and assume that the time-interval is [0, k], keeping in mind that, although not visible, all other components are present in f . We thus seek to determine the coefficients {ξ n } q n=1 with ξ 0 given, such that for m = 1, . . . , q we have n (τ (t))λ We explicitly compute the inverseĀ [ [q] mn a n0 = −1.
The weight functions may be computed analytically for small q, and for general q they are easy to compute numerically.
A.2. The mdG(q) method. We now make the same ansatz as for the continuous method, where the difference is that we now have q + 1 degrees of freedom on every interval, since we no longer have the continuity requirement for the trial functions. We make the assumption that the nodal points for the nodal basis functions are chosen so that for m = 0, . . . , q ij , where we use ξ − ij0 to denote ξ i,j−1,qi,j−1 , i.e., the value at the right end-point of the previous interval. To simplify the notation, we drop the subindices ij again and rewrite to [0, k]. We thus seek to determine the coefficients {ξ n } q n=0 such that for m = 0, . . . , q we have  mn ). Then, switching back to the full notation, we have (A.20) [q] mn λ [q] n , m = 0, . . . , q.
As for the continuous method, this may be somewhat simplified. Lemma A.2. For the mdG(q) method, we have q n=0ā [q] mn λ [q] n (0) = 1.
Proof. As in the proof for the mcG(q) method, assume that the interval is [0, 1]. Since the value of the expression is independent of f we can take f = 0. We thus want to prove that if f = 0, then the solution U is constant. By the Galerkin orthogonality, we have The mdG(q) method thus reads as follows: For every local interval I ij , find {ξ ijn } qij n=0 , such that for m = 0, . . . , q ij we have for certain weight functions {w [q] n } q n=0 ⊂ P q (0, 1).