# Introduction to Mathematical Physics/Some mathematical problems and their solution/Evolution, numerical time integration

## Introduction

There exist very good books (see ([#References|references])) about numerical integration of PDE evolution problem. Note that all methods (finite difference, finite element, spectral methods) contains as a final step the time integration of a ODE system. In this section this time integration is treated. Problem to be solved is the following

probmathevovec

Problem: (Cauchy problem) Find functions $u_k(t)$ obeying the ODE system

$\frac{d u_k}{dt}=F_k(u_h,t)$

where $u_h(x,t=0)$ is known.

We will not present here details about stability and precision calculation of the numerical integration schemes, however the reader should always keep in mind this crucial problem. Generally speaking, knowledge of mathematical properties of solutions should always be used to verify numerical results. For instance, if the solution of problem probmathevovec is known to be bounded, and that the solution obtained numerically diverges, then numerical solution will be obviously considered as bad. This problem of stability is the first that the numerician meet. However, more refined considerations have to be considered. Integrals of movement are a classical way to check accuracy of solutions. For hamiltonian systems for instance, energy conservation should be checked at regular time intervals.

## Euler method

Euler method consists in approximating the time derivative

$\frac{du_k}{dt}$

by

$\frac{1}{\Delta t}(u_{k,t+1}-u_{k,t}).$

However, the way the right-hand term is evaluated is very important for the stability of the integration scheme ([#References|references]). So, explicit scheme:

$\frac{1}{\Delta t}(u_{k,t+1}-u_{k,t})=F_i(u_t)$

is unstable if $F=\frac{\partial u}{\partial x}$. On another hand, implicit scheme:

eqimpli

$\frac{1}{\Delta t}(u_{k,t+1}-u_{k,t})=F_i(u_{t+1})$

is stable. This last scheme is called "implicit" because equation eqimpli have to be solved at each time step. If $F$ is linear, this implies to solve a linear system of equations (see ([#References|references])). One can show that the so-called Cranck-Nicholson scheme:

$\frac{1}{\Delta t}(u_{k,t+1}-u_{k,t}) = \frac{1}{2}[F_i(u_{t+1})+F_i(u_{t})]$

is second order in time. Another interesting scheme is the leap-frog scheme:

$\frac{1}{2\Delta t}(u_{k,t+1}-u_{k,t-1}) =F_i(u_{t})$

Example:

Let us illustrate Euler method for solving Van der Pol equation:

$\ddot{x}-\epsilon(1-x^2)\dot{x}+x=0$

which can also be written:

$\begin{matrix} \frac{dx}{dt}&=&y\\ \frac{dy}{dt}&=&\epsilon (1-x^2)y-x \end{matrix}$
File:Pol03i
$\epsilon=0.3$, initial conditions inside the attractor.} LABEL figpol03i
File:Pol03e
$\epsilon=0.3$, initial conditions out of the attractor.} LABEL figpol03e

Example: Let us consider the reaction diffusion system:

$\frac{\partial u}{\partial t}=\gamma(a-bu+\frac{u^2}{\nu})+\nabla^2u$

$\frac{\partial v}{\partial t}=\gamma(u^2-v)+d\nabla^2v$

Figure figreacd shows results of integration of this system using Euler method (space treatment is achieved using simple finite differences, periodic boundary conditions are used). The lattice has $100 \times 50$ sites.

\begin{figure} \begin{tabular}[t]{ccc}

\epsffile{reacdt1ga5e_04ite2} \epsffile{reacdt1ga5e_04ite4}

\epsffile{reacdt1ga5e_04ite6} \epsffile{reacdt1ga5e_04ite8}

\epsffile{reacdt1ga5e_04ite10} \epsffile{reacdt1ga5e_04ite12}

\epsffile{reacdt1ga5e_04ite14} \epsffile{reacdt1ga5e_04ite16}

\epsffile{reacdt1ga5e_04ite18} \epsffile{reacdt1ga5e_04ite20} \end{tabular}

Example:

Nonlinear Klein-Gordon equation

$u_{xx}-u_{tt}=F(u)$

where $F$ is (Sine-Gordon equation)

$F(u)=\sin u$

or (phi four)

$F(u)=-u+u^3$

can be integrated using a leap-frog scheme ([#References

## Runge-Kutta method

Runge-Kutta formula \index{Runge-Kutta} at fourth order gives $u_{i+1}$ as a function of $y_i$:

$u_{i+1}=u_i+\frac{\Delta t}{6}[f(u_i,t_i)+2f(uu_{i+\frac{1}{2}},t_{i+\frac{1}{2}})+ 2f(uuu_{i+\frac{1}{2}},t_{i+\frac{1}{2}})+ f(uu_{i+1},t_{i+1})]$

where

$uu_{i+\frac{1}{2}}=u_i+\frac{\Delta t}{2}f(u_i,t_i)$

$uuu_{i+\frac{1}{2}}=u_i+\frac{\Delta t}{2}f(uu_{i+\frac{1}{2}},t_{i+\frac{1}{2}})$

$uu_{i+1}=u_i+{\Delta t}f(uuu_{i+\frac{1}{2}},t_{i+\frac{1}{2}})$

Problem:

Find $u(t)$ solution of:

$\frac{d u}{dt}=F(u)$

with

$u(0)=u_0$

$u_{i+1}=y_i+\frac{\Delta t}{2}(3F_i-F_{i-1})$

or open Adams method order three:

$u_{i+1}=y_i+\frac{\Delta t}{12}(23F_i-16F_{i-1}+5F_{i-2})$

Open Adams method are multisteps: they can not initiate an integration process. Examples of closed Adams method, are at order two:

$u_{i+1}=y_i+\frac{\Delta t}{2}(F_{i+1}+F_{i})$

and order three:

$u_{i+1}=y_i+\frac{\Delta t}{12}(5F_{i+1}-8F_{i}-F_{i-1})$

Closed Adams methods are also multisteps. They are more accurate but are implicit: they imply the solving of equations systems at each time step. Predictor--corrector method provides a practical way to use closed Adams formulas.

## Predictor--corrector method

Predictor--corrector methods\index{predictor-corrector} are powerful methods very useful in the case where right--hand side is complex. An first estimation $\tilde u_{i+1}$ is done (predictor) using open Adams formula, for instance open Adams formula order three:

$\tilde u_{i+1}=u_i+\frac{\Delta t}{12}[23F(u_i)-16F(u_{i-1})+5F(u_{i-2})]$

Then, $u_{i+1}$ is re-evaluated using closed Adams formula where $\tilde u_{i+1}$ computed previously is used to avoid implicity:

$u_{i+1}=y_i+\frac{\Delta t}{12}[5F(\tilde u_{i+1})-8F(u_{i})-F(u_{i-1})]$

Remark:

Evaluation of $F$ in the second step is stabilizing. Problem

$\frac{du}{dt}=Lu+N(u)$

where $L$ is a linear operator and $N$ a nonlinear operator will be advantageously integrated using following predictor--corrector scheme:

$\tilde u_{i+1}=\exp(L.\Delta t)u_i+\Delta t N(u_i)$

$u_{i+1}=\exp(L.\Delta t)u_i+\Delta t N(\frac{\tilde u_{i+1}+u_i}{2})$

## Symplectic transformation case

Consider\index{symplectic} the following particular case of evolution problem:

Problem:

Find solution $[q(t),p(t)]$

$\frac{dq}{dt}=f(q,p)$

$\frac{dp}{dt}=g(q,p)$

where $[q(t=0),p(t=0)]$ is known (Cauchy problem) and where functions $f$ and $g$ come from a hamiltonian $H(q,p)$:

$f=\frac{\partial H}{\partial p}$

$g=-\frac{\partial H}{\partial q}$

Dynamics in this case preserves volume element $dq.dp$. Let us define symplectic \index{symplectic}transformation:

Definition:

Transformation $T$ defined by:

$\begin{matrix} q_{n+1}&=&f(q_n,p_n)\\ p_{n+1}&=&g(q_n,p_n) \end{matrix}$

is symplectic if its Jacobian has a determinant equal to one.

Remark:

eqsymple1

$\frac{dq}{dt}=V(p)$

eqsymple2

$\frac{dp}{dt}=F(q)$

consists in integrating using the following scheme:

$\begin{matrix} q_{n+1}&=&q_n+V(p_n)\delta t\\ p_{n+1}&=&p_n+F(q_n)\delta t \end{matrix}$

and is not symplectic as it can be easily checked.

Let us giev some examples of symplectic integrating schemes for system of equation eqsymple1 and eqsymple2.

Example: Verlet method:

$\begin{matrix} q_{n+1}&=&q_n+V(p_n)\delta t\\ p_{n+1}&=&p_n+F(q_{n+1})\delta t \end{matrix}$

is symplectic (order 1).

Example:

Staggered leap--frog scheme

$\begin{matrix} \frac{1}{\delta t}(p_{n+1}-p_{n})&=&F(q_{n+\frac{1}{2}})\\ \frac{1}{\delta t}(q_{n+1+\frac{1}{2}}-q_{n+\frac{1}{2}})&=&V(P_{n+1}) \end{matrix}$

is symplectic (order 2).

Figures codefigure1, codefigure2, and codefigure3 represent trajectories computed by three methods: Euler, Verlet, and staggered leap--frog method for (linear) system describing harmonic oscillator:

$\left\{ \begin{array}{lcr} \dot x&=&-y\\ \dot y&=&x \end{array}\right.$

\begin{figure}[tbh]

EulerOscillo001
| center | frame |\it Euler method with time-step $h=0.01$ for harmonic oscillator.}


LABEL codefigure1 ]]

\begin{figure}[tbh]

VerletOscillo1
| center | frame |\it Verlet method with time step $h=0.9$ for harmonic
oscillator. Despite a larger time-step than time step used to plot previous
figure, energy is well conserved.}


LABEL codefigure2 ]]

\begin{figure}[tbh]

VerletRecOscillo09
| center | frame |\it Staggered leap--frog method: energy is still well conserved, but
space phase is less deformed than in previous case.}


LABEL codefigure3 ]]