# 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 $u(t)$ such that

• $\dot u(t)+Au(t) = F(t)\quad\mbox{for }t>0,\quad u(0)=u^0$,

where $A$ is a constant and $F(t)$ a given functions, which has the standard form

• $\dot u(t)=f(t,u(t))$ with $f(t,u)=F(t)-Au$.

Compute an approximate solution $U(t)$ by time stepping according to the Trapezoidal Method (here the same as Midpoint Euler ): Find $U^{n}=U(nk)$ such that

• $U^{n+1}+\frac{k}{2}(AU^{n+1}+AU^n)=U^{n}+\int_{I_n}F(t)\, dt, \quad \mbox{for }n=0,1,2,...,$

with $I_n=(nk,(n+1)k)$ and $U^0\approx u^0$, assuming the integral of $F(t)$ can be evaluated analytically.

We here think of $U(t)$ as a piecewise linear continuous function taking on the computed values $U^n$ at the discrete time levels $nk$ with time step $k$. In this case the Trapezoidal Method and Midpoint Euler method coincide, since

• $\frac{1}{2}(U^{n+1}+U^n)=U(nk+\frac{k}{2})$,
where $nk + \frac{k}{2}$ is the midpoint between $nk$ and $(n+1)k$. We also have
• $\int_{I_n}\dot U\, dt = U^{n+1}-U^{n},\quad \int_{I_n}AU\, dt = \frac{k}{2}(AU^{n+1}+AU^n)$,

which shows that Trapezoidal/Midpoint Euler satisfies:

• $\int_{I_n}(\dot U+AU -F)\, dt =0\quad\mbox{for }n=0,1,2,...$

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

• $R(U)\equiv \dot U +AU - F$,
vanishes.

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

1. constant non-negative: $A\ge 0$,
2. constant imaginary: $A = i$,
3. constant negative: $A <0$,
4. oscillating positive-negative: $A(t)=\sin(t)$,

where we also added a basic case with $A(t)$ depending on $t$.

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

1. $A$ positive semi-definite (symmetric): diffusion problems
2. $A$ anti-symmetric: wave problems
3. $A$ negative definite (symmetric) or general matrix: inverted pendulum…
4. $A$ oscillating with zero mean: turbulence…

For a general non-linear system $\dot u+f(u)=0$, the matrix $A$ then corresponds to the Jacobian $f^\prime (u(t))$ as concerns stability.

# Error Analysis

We shall now estimate the error $\vert u(T)-U(T)\vert$ at a given time $T=Nk$ in terms of the time step $k$ and relevant quantities to be defined, where we assume the $u(t)$ 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 $\phi (t)$ of the following dual problem:

•  $-\dot \phi +A\phi =0 \mbox{ for } T> t \geq 0$, $\phi(T)=\pm 1= \mbox{sign of } e(T)$,             (1)

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

• $0=\int_{0}^{T} e\,(-\dot \phi +A\phi )\, dt$ ,

and integrate by parts to get the following error representation (since $\vert e(T)\vert =e(T)\phi (T)$):

• $\vert e(T)\vert=\int_{0}^{T}(\dot e +Ae)\phi\, dt+e(0)\phi (0)$,

where we allow $U(0)$ to be different from $u(0)$, corresponding to an error $e(0)$ in the initial value $u(0)$.

Since $u$ solves the differential equation, that is $\dot u+Au=F$, we have

• $\dot e+Ae=\dot u+Au-F-(\dot U+AU-F)=F-\dot U-AU=-R(U)$,

and thus we obtain the following representation of the error $\vert e(T)\vert$ in terms of the residual $R(U)=\dot U+AU-F$ and the dual solution $\phi$:

• $\vert e(T)\vert =-\int_{0}^{T}R(U)\phi\, dt +e(0)\phi (0)$.

Now Midpoint Euler states that

• $\int_{I_n}R(U)\,dt=0\quad\mbox{for }n=0,1,2,...,$

which allows us to rewrite the error representation as

• $\vert e(T)\vert =-\int_{0}^{T}R(U)(\phi-\bar\phi )\, dt+e(0)\phi (0)$,

where $\bar\phi$ is the mean-value of $\phi$ over each time interval $I_n$, that is

• $\bar\phi (t)=\frac{1}{k}\int_{I_n}\phi (s)\, ds\quad\mbox{for }t\in I_n$.

We shall now use the fact that

• $\int_{I_n}\vert\phi-\bar\phi\vert\, dt\le k\int_{I_n}\vert\dot \phi \vert\, dt$ ,

which follows by integration from the facts that

• $\phi (t)-\bar\phi (t)=\frac{1}{k}\int_{I_n}(\phi (t)-\phi (s))\, ds$,

and

• $\vert \phi (t)-\phi (s)\vert\le \int_{s}^t\vert \dot\phi (\sigma ) \vert \, d\sigma \le\int_{I_n}\vert \dot\phi (\sigma ) \vert \, d\sigma\quad\mbox{for }s, t\in I_n$.

We conclude that

• $\vert e(T)\vert\le \sum_{n=0}^{N-1}R_n\int_{I_n} \vert\phi-\bar\phi \vert dt +\vert e(0)\vert\vert\phi (0)\vert$
• $\le \sum_{n=0}^{N-1}k\, R_n\,\int_{I_n} \vert\dot \phi \vert dt +\vert e(0)\vert\vert\phi (0)\vert$ ,

where

• $R_n(U)=\max_{t\in I_n}\vert R(U(t))\vert$.

Bringing out the max of $k_nR_n$ over $n$, we get

• $\vert e(T)\vert\le\max_{0 \leq n \leq N-1}kR_n\int_0^T\vert\dot\phi\vert dt+\vert e(0)\vert\vert \phi(0)\vert$.

Defining the stability factors $S_c(T)$ and $S_d(T)$ by

• $S_c(T)=\int_0^T\vert\dot\phi (s)\vert\, ds,\quad S_d(T)=\vert\phi (0)\vert$,

we get the following a posteriori error estimate:

Theorem: The approximate solution $U(t)$ of the initial value problem $\dot u(t)+Au(t) = F(t)\quad\mbox{for }t>0,\quad u(0)=u^0$, computed by Midpoint Euler with time step $k$ over intervals $I_n$, satisfies for $T>0$

• $\vert u(T)-U(T)\vert\le S_c(T)\max_{n}kR_n(U) +S_d(T)\vert u(0)-U(0)\vert$,

where $u(t)$ is the exact solution computed with vanishingly small time step, $R_n(U)=\max_{I_n}\vert\dot U+AU-F\vert$ measures the residual over $I_n$, and $S_c(T)$ and $S_c(T)$ are stability factors related to time-discretization and initial data.

The stability factors $S_c(T)$ and $S_d(T)$ 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 $S_c(T)$ and $S_d(T)$ depending on the nature of $A$,in particular the sign of $A$, with possibly vastly different stability factors.

We notice that the solution $\phi (t)$ of  the dual problem in the case $A$ is constant, is given by the explicit formula

• $\phi (t)=\pm \exp(-A(T-t))$.

We see that if $A >0$, then the solution $\phi (t)$ decays as $t$ decreases from $T$,and the case $A > 0$ is thus the “stable case”. If $A<0$ then the exponential.factor $\exp(-AT)$ 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 $S_c(T)$ and $S_d(T)$ satisfy:

• $S_d(T) \le \exp(\vert A\vert T),\quad S_c(T)\le \exp(\vert A\vert T)$ if $A<0$ ,
• $S_d(T) \le 1,\quad S_c(T)\le 1$ if $A\ge 0$ ,
• $S_d(T) \le 1,\quad S_c(T)\le T$ if $A=i$ ,
• $S_d(T) \le \exp(1),\quad S_c(T)\le exp(1)T$ if $A=\sin(t)$ .

Proof: Changing variables $T-t\rightarrow t$, we can write the dual equation as the forward-in-time problem $\dot\phi =-A\phi$ for $t>0$, $\phi (0)= 1$ with solution $\exp(-At)$, if $A$ is constant. We note that if $A>0$,

• $\int_0^T\vert \dot\phi (t)\vert\, dt =\int_0^TA\exp(-At)\, dt = -\int_0^T\frac{d}{dt}\exp(-At)\, dt=1-\exp(-AT)<1$.

Further, if $A<0$, then

• $\int_0^T\vert \dot\phi (t)\vert\, dt =\int_0^T\frac{d}{dt}\exp(-At)\, dt=\exp(-AT)-1\approx \exp(\vert A\vert T)$

If $A=i$, then

• $\int_0^T\vert \dot\phi (t)\vert\, dt =\int_0^T\vert\frac{d}{dt}\exp(-it)\vert\, dt=\int_0^T1\, dt=T$.

Finally, if $A=\sin(t)$ then $\phi (t)=\exp(\cos(t))$, and so

• $\vert\phi(T)\vert\le\exp (1),\quad \int_0^T\vert\dot\phi (t)\vert\, dt\le\exp(1)T$.

$\blacksquare$

The size of the stability factors indicate the degree of stability of the solution $u(t)$ being computed. If the stability factors are large, the residuals $R(U(t))$ and $e(0)$ 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 $U(t)$. There is a corresponding \emph{a priori error estimate} with $R(U)$ replaced by $R(u_k)$ where $u_k$ is the piecewise linear interpolant of the exact solution $u(t)$ taking on the same values at the discrete time levels $nk$. In this case the stability factors measure the stability of a corresponding \emph{discrete dual problem}.

How big is then $R(u_k)$? Well, with piecewise linear interpolation, we have $\vert\dot u -\dot u_k\vert \approx k\vert \ddot u\vert$, and thus the a priori estimate takes the form

• $\vert u(T)-U(T)\vert\le S_c(T)C(u)k^2 +S_d(T)\vert e(0)\vert$,

where $C(u)=max_t\vert \ddot u(t)\vert$. In short, Midpoint Euler is second-order accurate with error proportional to $k^2$. Backward Euler and Forward Euler are first order accurate witherror proportional to $k$.

# Generalization

The above error analysis extends to a general IVP $\dot u(t)=f(u(t))$ for $t>0$, $u(0)=u^0$, as shown in Chapter (\ref{chapterFEMIVP}).

• Show that the above posteriori error estimate directly extends to variable time steps $k_n$ with $k_nR_n$ replacing $kR_n$.