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}).

To Think About

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: