I believe in the fundamental Truth of all the great religions of the world. I believe that they are all God given. I came to the conclusion long ago… that all religions were true and also that all had some error in them. (Mahatma Gandhi)
The fairest thing we can experience is the mysterious. It is the fundamental emotion which stands at the cradle of true art and true science. He who know it not and can no longer wonder, no longer feel amazement, is as good as dead, a snuffed-out can. (Einstein)
Most of the fundamental ideas of science are essentially simple, and may, as a rule, be expressed in a language comprehensible to everyone. (Einstein)

Time Stepping

The Fundamental Theorem concerns time-stepping of the basic IVP (initial value problem):

$Du(t)=f(t)\quad\mbox{for }t>0,\quad u(0)=u_0$ ,

where the Lipschitz continuous function $f(t)$ does not depend on the unknown $u$ , only on the (independent) time variable t, $Du = \frac{du}{dt}$ and $u_0$ is a given initial value.

We now extend to allow the function $f(\cdot )$ to depend also on $u$ . Assuming for simplicity no explicit dependence on t, we thus consider the IVP:

$Du(t)=f(u(t))\quad\mbox{for }t>0,\quad u(0)=u_0$ ,

where $f(u)$ is a given Lipschitz continuous function of $u$ with Lipschitz constant $L$ .

The basic question is if the solution can be computed to arbitrary precision by time-stepping $du = f(u)*dt$ ?

The basic case is $f(u)=u$ and $u_0=1$ , that is the IVP:

$Du(t)=u(t)\quad\mbox{for }t>0,\quad u(0)=1$ ,

with $L=1$ and solution denoted by $u(t)=\exp(t)$ referred to as the exponential function being computed by Forward Euler:

$u(t+dt) = u(t) + dt*u(t) = (1+dt)*u(t)$ ,
$u(t+n*dt) = (1+dt)^n$ ,
$\exp(t)\approx (1+\frac{t}{n})^n$ ,

with time step $dt = t/n$ .

We estimate the effect of dividing the time-step by a factor 2, using that $(t+dt)^n-t^n\approx nt^{n-1}dt$ (because $\frac{d}{dt}t^n=nt^{n-1}$ ):

$(1+\frac{t}{2n})^{2n}-(1+\frac{t}{n})^{n}=((1+\frac{t}{2n})(1+\frac{t}{2n}))^n-(1+\frac{t}{n})^n$
$=((1+\frac{t}{n}+\frac{t^2}{4n^2})^n-(1+\frac{t}{n})^n\approx n(1+\frac{t}{n})^{n-1}\frac{t^2}{4n^2}$
$\approx \frac{t}{n}\exp(t)\frac{t}{4}$ .

We see that the difference is proportional to the time step $\frac{t}{n}$ . As in the proof of the Fundamental Theorem of Calculus, we conclude that $(1+\frac{t}{n})^n$ determines $\exp(t)$ with a precision proportional to the time step with a multiplicative factor $\approx \exp(t)\frac{t}{4}\sim \exp(t)$ .

Time Stepping with General f(u)

The above proof extends to an arbitrary Lipschitz continuous function $f(u)$ with the difference that the time-stepping error in computing $u(t)$ , now is proportional to the time step with a multiplicative factor $\exp(L*t)$ , where $L$ is the Lipschitz constant of $f(u)$ . This extends to systems with $f:R^d\rightarrow R^d$ with d>1.

It is natural to refer to this result as a Generalized Fundamental Theorem: The Fundamental Theorem concerns $Du=f(t)$ and the Generalized Fundamental Theorem concerns $Du=f(u(t))$ . Calculus in a nutshell!

The Two Parts of the Proof

A proof of the Generalized Fundamental Theorem can be performed by combining the following two parts:

Part 1: Estimate the difference $u((n+1)k)-\tilde u((n+1)k)$ by taking one (Forward Euler) time step of length $k$ (here $dt = k$ ) and two time steps of length $k/2$ , from the same initial value $u(nk)$ :

$\vert u((n+1)k)-\tilde u((n+1)k)\vert$
$=\vert u(nk)+kf(u(nk))-(u(nk)+\frac{k}{2}f(u(nk))+\frac{k}{2}f(u(nk)+\frac{k}{2}f(u(nk)))\vert$
$=\frac{k}{2}\vert f(u(nk))-f(u(nk)+\frac{k}{2}f(u(nk)))\vert$
$\le\frac{k}{2}L\frac{k}{2}\vert f(u(nk))\vert$ .

Illustration of Part 1 in the proof of the Generalized Fundamental Theorem

Part 2: Estimate the difference after one time step from different initial conditions $u(nk)-\tilde u(nk)$ :

$\vert u((n+1)k)-\tilde u((n+1)k)\vert =\vert u(nk)+kf(u(nk))-\tilde u(nk)+kf(\tilde u(nk))\vert$
$\le (1+kL)\vert u(nk)-\tilde u(nk)\vert$ .

Illustration of Part 2 of the proof of the Generalized Fundamental Theorem.

Combining (see hint below) Part 1 and Part 2, we obtain a final error proportional to the time step $k$ with a multiplicative factor $\exp(Lt)$ , which we refer to as a stability factor.

To see the connection with the basic case $f(u)=u$ , think of estimating a general function $f(u)$ , assuming $f(0)=0$ for simplicity, by $f(u)\approx f^\prime (0)*u$ , which suggests that the factor $\exp(t)$ for $f(u)=u$ should be replaced by $\exp(L*t)$ for a general $f(u)$ (because $\vert f^\prime (0)\vert \le L$ ).

Generalization to a vector valued function $f:R^d\rightarrow R^d$ with $d>1$ is direct, and we have thus presented the essentials of a proof (with a hint to completion below) of the following main result of Calculus:

Generalized Fundamental Theorem of Calculus: The solution u(t) of the IVP $\dot u(t)=f(u(t))$ for $0<t\le T$ with u(0) given, where $f:R^d\rightarrow R^d$ is Lipschitz continuous with Lipschitz constant $L$ , is uniquely computable by Forward Euler time stepping with a precision proportional to the time step times a stability factor of size exp(LT).

The proof is similar for the alternative methods (with $u^n=u(nk)$):

$u((n+1)k)=u(nk)+kf(u(nk))$ Forward Euler
$u^{n+1}=u^n+kf(u^{n+1})$ Backward Euler
$u^{n+1}=u^n+kf(\frac{u^n+u^{n+1}}{2})$ Midpoint Euler
$u^{n+1}=u^n+\frac{k}{2}(f(u^n)+f(u^{n+1}))$ Trapezoidal Method.

Note that if $f(u)$ is linear in $u$ , then Midpoint Euler and the Trapezoidal Method coincide.

A Posteriori Error Control

For a more precise error control, based on computed solutions, see

MST Chap 82: Time Stepping Error Analysis,
MST Chap 161: Time Stepping by FEM.

The Illusion of an Exponential Bound

If $L=10$ and $T=30$ , which looks pretty harmeless, then $\exp(LT)=\exp(300) >>10^{100}$ = googol, an incredibly large number, much larger than the number of atoms in the Universe.

A matching time step of $10^{-100}$ is beyond all rationale and thus computation of a solution of an IVP with moderate Lipschitz constant over a moderatley long time interval may be impossible. An example is the Lorenz system with

$f(u)=(-10 u_1+10 u_2, 28u_1-u_2-u_1u_3,-\frac{8}{3}u_3+u_1u_2)$ ,

for which computation on an interval of length $T$ requires computation with about $T/2$ digits, see:

MST Chap 256: The Lorenz System and the Essence of Chaos
Long-Time Computability of the Lorenz System

Stiff IVPs

There is a class of IVPs with large or very large Lipschitz constants, which are computable on long time intervals, because the function $f(u)$ has a decay property (negative derivative) causing errors to decay rather than grow exponentially. Such problems are called stiff problems and may require implicit time stepping to avoid severe time step restrcitions in explicit methods. See \hyperref[sectionstiffproblems]{Stiff Problems}.

lorenz — A Lorenz system solution trajectory.

Wave Equations

IVPs with wavelike solutions, like the system

$\dot u_1=u_2\quad \dot u_2=-u_1$

with solutions being linear combinations of $\sin(t)$ and $\cos(t)$ , have formally Lipschitz constants of size 1, can be integrated with error growth $\sim t$ , instead of $\sim\exp(t)$ by the above (crude) estimate, if a proper time-stepping method like Trapezoidal method is used. This is due to error cancellation in wave motion.

For more complex wave problems, or problems with more or less periodic solutions, the stability factorcan have a polynomial growth in time $t$ , e.g. quadratic for simple planetary systems.

Summary: Time Stepping of IVP

The precision in time stepping the solution $u(t)$ of an IVP $\dot u=f(u)$ for $0<t\le T$, with a first order method such as Forward/Backward Euler with time step $k$ , can be estimated by $S(T)k$ , where $S(t)$ acts as stability factor measuring error propagation and accumulation of size

$S(T)\sim\exp(LT)$ (general), where $L$ is the Lipschitz constant of $f$,
$S(T)\sim 1$ (stiff: diffusion problems)
$S(T)\sim T$ (wave problems), $S(T)\sim T^2$ (planetary system).

Preparing for a More Precise Analysis

As a preparation for the more precise error analysis in

MST Chap 82: Time Stepping Error Analysis,
MST Chap 161: Time Stepping by FEM,

we consider two solutions $u(t)$ and $\tilde u(t)$ computed with the same time step $k$ but different initial data $u^0$ and $\tilde u^0$ . Subtracting the update formulas we have formally for the difference $e=u-\tilde u$ :

$\dot e(t)\approx f^\prime (u(t))e(t)\quad\mbox{for }t>0,\quad e(0)=e^0\equiv u^0-\tilde u^0$

showng that an initial error is propagated as a solution to a linearized IVP with coefficient $f^\prime (u(t))$ depending on a computed solution $u(t)$ . We shall see that by solving the linearized problem (or rather a closely related dual linearized problem), the stability factors $S(t)$ measuring error growth can be computed and the precision on the computation of $u(t)$ can be assessed.

The linearized problem (or its dual) thus gives the key to unlock time stepping precision.

Note that with $f(u)=u$ , the linearized problem reads $\dot e=e$ with solution $e(t)=\exp(t)e^0$ , showing exponential error growth, as expected.

Completion of the Proof

To complete the proof of the Generalized Fundmental Theorem we are to sum up the error contributions from each subinterval of length $k$ , which according to Part 1 and Part 2 amounts to

$\sum_{n=1}^N(1+kL)^nLk^2M\approx kML\sum_{n=0}^N\exp(Lnk)k$
$\approx kML\int_0^T\exp(Ls)\, ds\approx kM\exp(LT)$ ,

where $Nk=T$ and $M\ge\max_{u}\vert f(u)\vert$. Can you explain what is going on here? If not, take a look at:

Hint to Completion of the Proof

We can think of comparing computations with $k$ and $\frac{k}{2}$ with corresponding solutions $u(t)$ and $\tilde u(t)$ in two ways depending on how we choose initial values on each time interval $(nk,(n+1)k)$ :

Compute $u(t)$ and $\tilde u(t)$ independently with timestep $k$ and $\frac{k}{2}$ .
Assume that $\tilde u(nk)=u(nk)$ at the beginning of time step $n$ and account for the effect at final time of the difference $\tilde u((n+1)k)-u((n+1)k)$ at the end of time step $n$ .

1. is the most direct from computational point of view and a corresponding proof is given in \hyperref[chaptergenivp]{General Case}.

Here we consider 2. because the proof is (maybe) simpler: The error from the first time step is bounded by $Lk^2M$ assuming no error in initial data, and is propagated with a factor bounded by $(1+kL)$ for each time step, thus with a factor $(1+kL)^N$ after $N$ steps. Similarly, the error from the second time step is bounded by $Lk^2M$ , again assuming no error in the corresponding initial value, and is propagated with a factor $(1+kL)^{N-1}$ , et cet. Summing we obtain a bound of the total error after $N$ steps.

Uniqueness of Solution

To prove that the solution of the $\dot u=f(u)$ with $u(0)=u^0$ is unique, we assume $v(t)$ is a possibly different solution also satisfying $\dot v(t)=f(v(t))$ for $t>0$ and $v(0)=u^0$ . Subtraction gives for the difference $w=u-v$

$\dot w=f(u)-f(v)$

and thus taking the scalar product with $w$ and using Cauchy’s inequality, we get

$\frac{d}{dt}\frac{1}{2}\vert w\vert^2 =\frac{d}{dt}\frac{1}{2}w\cdot w=(f(u)-f(v))\cdot w\le L\vert w\vert\vert w\vert$

and thus for $W=\vert w\vert^2$

$\dot W\le 2LW$ ,

which shows that (why?)

$W(t)\le W(0)\exp(2Lt)\quad\mbox{for }t>0$ .

But $W(0)=\vert u(0)-v(0)\vert =0$ and thus $W(t)=0$ and $u(t)=v(t)$ for $t>0$ and uniqueness follows. But to be scientifically honest, size of the exponential factor $\exp(Lt)$ is crucial. If $L=10$ and $t=30$ , which does not look too frightening, then $\exp(Lt)=\exp(300)>10^{100}=$ googol, which means that the argument the $\exp(300)W(0)$ is small (zero) if $W(0)$ is small (zero), is questionable, very questionable, right?

How to Prove that exp(t+s)=exp(t)exp(s)?

To prove the basic law of the exponential $\exp(t+s)=\exp(t)\exp(s)$ , note that the function $u(s)=\exp(t+s)$ satisfies $\frac{du}{ds}=u$ for $s>0$ and $u(0)=\exp(t)$ . But the function $v(s)=\exp(t)\exp(s)$ also satisfies $\frac{dv}{ds}=v$ for $s>0$ and $v(0)=\exp(t)$ and by uniqueness $u(s)=v(s)$ , that is $\exp(t+s)=\exp(t)\exp(s)$ .