# Time Stepping Error Analysis

*Sapiens nihil affirmat quod non probat.**After experience had taught me that all the usual surroundings of social life are vain and futile; seeing that none of the objects of my fears contained in themselves anything either good or bad, except in so far as the mind is affected by them, I finally resolved to inquire whether there might be some real good having power to communicate itself, which would affect the mind singly, to the exclusion of all else: whether, in fact, there might be anything of which the discovery and attainment would enable me to enjoy continuous, supreme, and unending happiness.*(Spinoza)

# Trapezoidal Method/Midpoint Euler

Consider a linear scalar IVP of the form: Find such that

- ,

where is a constant and a given functions, which has the standard form

- with .

Compute an approximate solution by time stepping according to the Trapezoidal Method (here the same as Midpoint Euler ): Find such that

with and , assuming the integral of can be evaluated analytically.

We here think of as a piecewise linear continuous function taking on the computed values at the discrete time levels with time step . In this case the Trapezoidal Method and Midpoint Euler method coincide, since

- ,

- ,

which shows that Trapezoidal/Midpoint Euler satisfies:

In other words, the mean-value over each subinterval of the *residual*

- ,

We shall consider the following basic choices of $A$ with different stability characteristics connecting to \hyperref[sectiontimesteppingIVP]{Summary: Timestepping of IVP}:

- constant non-negative: ,
- constant imaginary: ,
- constant negative: ,
- oscillating positive-negative: ,

where we also added a basic case with depending on .

We shall below see that (\ref{scalinode1}) can also be interpreted as a system of differential equations with a square matrix with the following analogous stability characteristics:

- positive semi-definite (symmetric): diffusion problems
- anti-symmetric: wave problems
- negative definite (symmetric) or general matrix: inverted pendulum…
- oscillating with zero mean: turbulence…

For a general non-linear system , the matrix then corresponds to the Jacobian as concerns stability.

# Error Analysis

We shall now estimate the error at a given time in terms of the time step and relevant quantities to be defined, where we assume the is an exact solution satisfying (\ref{scalinode1}) to high precision (computed with a (vanishingly) small time step).

We shall do this by representing the error in terms of the solution of the following *dual problem*:

- , , (1)

where . We note that that (1) runs “backwards in time” starting at time , because the (initial) data is given at , and that (accordingly) the time derivative term $\dot \phi$ has a minus sign. We start from the identity

- ,

and integrate by parts to get the following error representation (since ):

- ,

where we allow to be different from , corresponding to an error in the initial value .

Since solves the differential equation, that is , we have

- ,

and thus we obtain the following representation of the error in terms of the residual and the dual solution $\phi$:

- .

Now Midpoint Euler states that

which allows us to rewrite the error representation as

- ,

where is the mean-value of over each time interval , that is

- .

We shall now use the fact that

- ,

which follows by integration from the facts that

- ,

and

- .

We conclude that

- ,

where

- .

Bringing out the max of over , we get

- .

Defining the *stability factors* and by

- ,

we get the following *a posteriori error estimate*:

**Theorem:** The approximate solution of the initial value problem , computed by Midpoint Euler with time step over intervals , satisfies for

- ,

where is the exact solution computed with vanishingly small time step, measures the residual over , and and are stability factors related to time-discretization and initial data.

The stability factors and measure the effects of the accumulation of error in theapproximation. To give the analysis a quantitative meaning, we have to give a quantitative bound ofthese factors. In general the stablity factors are computed by computing the solution of the dual problem.

In special cases the stability factors can be computed analytically, as we now show: The following lemma gives an estimate for and depending on the nature of ,in particular the sign of , with possibly vastly different stability factors.

We notice that the solution of the dual problem in the case is constant, is given by the explicit formula

- .

We see that if , then the solution decays as decreases from ,and the case is thus the “stable case”. If then the exponential.factor enters, and depending on the size of $A$ this case is “unstable”.

More precisely, we conclude directly from the explicit solution formula that

**Theorem**: The stability factors and satisfy:

- if ,
- if ,
- if ,
- if .

**Proof**: Changing variables , we can write the dual equation as the forward-in-time problem for $t>0$, with solution , if is constant. We note that if ,

- .

Further, if , then

If , then

- .

Finally, if then , and so

- .

The size of the stability factors indicate the degree of stability of the solution being computed. If the stability factors are large, the residuals and have to be made correspondingly smaller by choosing smaller time steps and the computational problem is more demanding.

**A Priori Error Estimate**

The a posteriori error estimate (\ref{apostivp}) estimates the error in terms of the computed solution . There is a corresponding \emph{a priori error estimate} with replaced by where is the piecewise linear interpolant of the exact solution taking on the same values at the discrete time levels . In this case the stability factors measure the stability of a corresponding \emph{discrete dual problem}.

How big is then ? Well, with piecewise linear interpolation, we have , and thus the a priori estimate takes the form

- ,

where . In short, Midpoint Euler is *second-order accurate *with error proportional to . Backward Euler and Forward Euler are *first order accurate* witherror proportional to .

# Generalization

The above error analysis extends to a general IVP for , , as shown in Chapter (\ref{chapterFEMIVP}).

# To Think About

- Show that the above posteriori error estimate directly extends to variable time steps with replacing .
- For a basic aspect of duality in error estimation, see \hyperref[chaptererrorduality]{Error Control by Duality}.

## Leave a Reply