previous up next ','..','$myPermit') ?> SYLLABUS  Previous: 1.2 Differential Equations  Up: 1.2 Differential Equations  Next: 1.2.2 Partial differential equations


1.2.1 Ordinary differential equations

Slide : [ ODE Euler - mid-point rule - Runge-Kutta || VIDEO 99) echo " modem - ISDN - LAN "; else echo "login";?> ]

Rather than repeating details of elementary numerical analysis, this section is only meant to review a few methods that will be used later, e.g. in the chapter dealing with the Monte Carlo method for partial differential equations (PDE).

Ordinary differential equations (ODE) are often encountered when dealing with initial value problems. In Newtonian mechanics, for example, the trajectory of a particle is described by an evolution in time

$\displaystyle \frac{d^2 \mathbf{X}}{dt^2}=\frac{F}{m} \quad\quad\Longleftrighta... ...ay}\right) = \left(\begin{array}{c}\mathbf{V} \\ \mathbf{F}/m\end{array}\right)$    

The example shows how higher order equations can be recast into a system of first order equations, with components of the form

$\displaystyle y^\prime = f(t,y), \quad\quad y(0)=y_0$ (1.2.1#eq.1)

Under very general assumptions, an initial condition $ y(0)=y_0$ having the same number of components as the unknown vector $ y$ yields exactly one solution to (1.2.1#eq.1), which can then be approximated with a computer.

Introducing a discretization $ y_n=y(t_n)$ with a finite number $ n$ of time steps $ t_n=n h +t_0$ , a straight forward manner to solve an ODE is to approximate the derivative with a finite difference quotient forward in time $ y^\prime=(y_{n+1}-y_{n})/h +\mathcal{O}(h^2)$ , leading to the Euler recursion formula

$\displaystyle y_{n+1} = y_n +h f(t_n,y_n), \quad\quad y_0=y(t_0)$ (1.2.1#eq.2)

Because all the quantities are known at time $ t_n$ , the scheme is said to be explicit. An implicit evaluation of the function $ f(t_{n+1},y_{n+1})$ is sometimes desirable in order to stabilize the propagation of numerical errors; this is however computationally expensive when the function cannot be inverted analytically.

More precision can be obtained by using a finite difference quotient symmetric in time $ y^\prime=(y_{n+1}-y_{n-1})/2h +\mathcal{O}(h^4)$ , which is obtained by summing two Taylor series forward and backward (exercise 1.02). This yields the so-called midpoint formula

$\displaystyle y_{n+1} = y_{n-1} +2h f(t_n,y_n),\quad\quad y_{-1}=y_0 -h f(t_0,y_0)$ (1.2.1#eq.3)

which is second order accurate, but requires a special initialization to generate the additional values that are needed from the past. Writing the mid-point rule as

$\displaystyle y_{n+1} = y_n +h f\left(t_n+\frac{h}{2}, y(t_n+\frac{h}{2})\right)$    

the initialization problem is cured and leads to the famous second order Runge-Kutta method
$\displaystyle k_1$ $\displaystyle =$ $\displaystyle hf(t_n,y_n)$  
$\displaystyle k_2$ $\displaystyle =$ $\displaystyle hf(t_n+\frac{h}{2}, y_n+\frac{1}{2}k_1)$  
$\displaystyle y_{n+1}$ $\displaystyle =$ $\displaystyle y_n+k_2$ (1.2.1#eq.4)

Such elementary and more sophisticated methods are commonly available in software packages such as OCTAVE and MATLAB; type help ode23 to find out more about how to use them. Methods for ordinary equations are usually robust; however, they become extremely inefficient when a problem is stiff, i.e. involves two very different scales by limiting the step size to the shorter scale even when this is not at all relevant (exercise 1.06).

SYLLABUS  Previous: 1.2 Differential Equations  Up: 1.2 Differential Equations  Next: 1.2.2 Partial differential equations