next up previous
Next: Seeing the flight path Up: The Art of Phugoid Previous: Existence of Solutions

Subsections

Numerical Methods

What does Maple do when you ask it to display a solution of a differential equation with DEplot? Obviously, it cannot solve the equations analytically, since solutions don't always exist in closed form. Instead, it approximates them numerically.

Euler's method

To get an idea of how this can be done, take a look again at the direction field for the glider. By looking at how the arrows point, you can almost guess how the solution must turn. From a given initial condition, move in the direction of the arrow to the next position. At that new point, step in the direction of the vector at that point, and so on. This is the idea behind the simplest numerical integration scheme, called Euler's method.

\begin{mfigure}\centerline{\psfig{figure=glider02.eps,height=1.5in}}\end{mfigure}


More precisely, given a differential equation and an initial condition

$\displaystyle {\frac{d\vec{x}}{dt}}$ = $\displaystyle \vec{F} $($\displaystyle \vec{x} $, t)                $\displaystyle \varphi$(t0) = $\displaystyle \vec{x}_{0}^{}$,

the goal is to construct a sequence of points $ \vec{x}_{0}^{}$,$ \vec{x}_{1}^{}$,$ \vec{x}_{2}^{}$,...,$ \vec{x}_{n}^{}$ which approximate the solution $ \varphi$ at times t0, t1, t2,..., tn. We take our initial point to be the initial condition, and choose some small number h > 0 called the stepsize. Then we set

$\displaystyle \vec{x}_{i+1}^{}$ = $\displaystyle \vec{x}_{i}^{}$ + h$\displaystyle \vec{F} $($\displaystyle \vec{x}_{i}^{}$, ti)                ti + 1 = ti + h

for 0$ \le$i$ \le$n - 1. This very simple scheme works best for small values of h; the error in the approximation for a fixed time interval is proportional to h. This means that to improve the accuracy of a numerical solution by a factor of 10, we need to do about 10 times more work.

Maple doesn't typically use Euler's method because the errors accumulate too quickly. To see this, we reproduce one of the earlier figures from §2 using Euler's method.

> 
  DEplot( phug, [theta(t), v(t)], t=0..5, 
      [seq([theta(0)=0, v(0)=1+0.2*i], i=0..10)],
      theta=-Pi..3*Pi, v=0.01..3, linecolor=black, method=classical[foreuler]);

\begin{mfigure}\centerline{ \psfig{figure=glider04.eps,height=2in}}\end{mfigure}

While the errors in the above figures are considerable, using a much smaller stepsize (say, stepsize=0.01) will essentially reproduce the earlier figure, with a lot more computation.


Numerical Solutions, Numerical Integration, and Runge-Kutta

To get some ideas about improving on Euler's method, let's first notice that Euler's method corresponds exactly to the left-hand Riemann sum when computing a numerical approximation to an integral. You may recall from calculus that the left-hand sum converges to the value of the integral, but you must take a large number of subdivisions to get any accuracy.


A more efficient method is the trapezoid rule, which is the average of the left-hand and right-hand sum. There is a numerical scheme for integrating ODEs which corresponds to this, known variously as the improved Euler method or the Heun formula. In the trapezoid rule, we compute the function at the left and right endpoints of the function and average the result; the analog of this would require us to average the value of the vector field at our current point and at the next point. Perhaps you see the circularity implicit in that reasoning: to know the next point, we would need to solve the differential equation. So, instead we fudge a bit: we get our approximation of the next point by using Euler's method. More specifically, using the notation of §4.1, given an approximation $ \vec{x}_{i}^{}$ at time ti, we find the next one using the scheme below.

ti + 1 = ti + h    
$\displaystyle \vec{k}_{\ell}^{}$ = $\displaystyle \vec{F} $($\displaystyle \vec{x}_{i}^{}$, ti)          slope at the left side of the interval
$\displaystyle \vec{k}_{r}^{}$ = $\displaystyle \vec{F} $($\displaystyle \vec{x}_{i}^{}$ + h$\displaystyle \vec{k}_{\ell}^{}$, ti + 1)          slope at the right side of the interval
$\displaystyle \vec{x}_{i+1}^{}$ = xi + h($\displaystyle \vec{k}_{\ell}^{}$ + $\displaystyle \vec{k}_{r}^{}$)/2    

To do the improved Euler method, we need to evaluate the vector field twice as often as the regular Euler method, but there is a big gain in accuracy. The error of the approximation for a fixed time interval is proportional to h2. Thus, if we decrease h by a factor of 10, we should expect to reduce the error 100-fold. You can get Maple to use the improved Euler method by specifying method=classical[heunform] as an argument to DEplot and related commands.


You might remember from your calculus course that Simpson's rule can get very accurate numerical approximations for integrals with a very few number of points. This is because Simpson's rule is a fourth-order method: the error in the approximation is proportional to h4. For Simpson's rule, we evaluate the function at the right endpoint, the left endpoint, and the middle, and then take a weighted average of those values.

The analog of Simpson's rule for differential equations is called the Runge-Kutta method3.2, developed by the German mathematicians C. Runge and W. Kutta at the end of the nineteenth century. It takes a linear combination of four slopes, as below:

ti + 1 = ti + h    
$\displaystyle \vec{k}_{\ell}^{}$ = $\displaystyle \vec{F} $($\displaystyle \vec{x}_{i}^{}$, ti)          slope at the left side of the interval
$\displaystyle \vec{k}_{m}^{}$ = $\displaystyle \vec{F} $($\displaystyle \vec{x}_{i}^{}$ + $\displaystyle {\frac{h}{2}}$$\displaystyle \vec{k}_{\ell}^{}$, ti + $\displaystyle {\frac{h}{2}}$)   slope at the middle of the interval
$\displaystyle \vec{k}_{c}^{}$ = $\displaystyle \vec{F} $($\displaystyle \vec{x}_{i}^{}$ + $\displaystyle {\frac{h}{2}}$$\displaystyle \vec{k}_{m}^{}$, ti + $\displaystyle {\frac{h}{2}}$)   corrected slope at the middle of the interval
$\displaystyle \vec{k}_{r}^{}$ = $\displaystyle \vec{F} $($\displaystyle \vec{x}_{i}^{}$ + h$\displaystyle \vec{k}_{c}^{}$, ti + 1)          slope at the right side of the interval
$\displaystyle \vec{x}_{i+1}^{}$ = xi + h($\displaystyle \vec{k}_{\ell}^{}$ + 2$\displaystyle \vec{k}_{m}^{}$ + 2$\displaystyle \vec{k}_{c}^{}$ + $\displaystyle \vec{k}_{r}^{}$)/6    

When the slope depends only on t and not on $ \vec{x} $, then $ \vec{k}_{m}^{}$ and $ \vec{k}_{c}^{}$ are the same, and the formula reduces exactly to that of Simpson's rule.

Runge-Kutta is probably the most commonly used method of numerical integration, because of its generally high accuracy. Like Simpson's rule, it is a fourth-order method, so the error is proportional to h4. Maple has several numerical methods for ODEs built in to it; see the help page on dsolve[numeric] for more information about them; the ones we have described are ``classical'' methods, and are described (along with others) on Maple's help page for dsolve[classical]. Unless asked to do otherwise, Maple's DEplot command uses the fourth-order Runge-Kutta method described above (Maple's name for this method is classical[rk4]), although the numeric option of dsolve defaults to a variation called fourth-fifth order Runge-Kutta-Fehlberg, which Maple refers to as rkf45.



Footnotes

... method3.2
In fact, there is a whole family of such methods, of various orders. But generally, if someone says ``Runge-Kutta'' without specifying the order, they mean this fourth-order Runge-Kutta.

next up previous
Next: Seeing the flight path Up: The Art of Phugoid Previous: Existence of Solutions

Translated from LaTeX by Scott Sutherland
2002-08-29