Numerical integration of ode

From HandWiki


Let Hepa img757.gif , yn (x)) be n functions of one variable x, with Hepa img758.gif and Hepa img759.gif the first and second derivatives. A first order, x-independent ODE has the form

Hepa img760.gif

A second order, x-dependent ODE has the form

Hepa img761.gif

In principle, these two forms are completely equivalent, the one is a special case of the other (if Hepa img762.gif and Hepa img763.gif , then Hepa img764.gif ; if Hepa img765.gif and Hepa img766.gif then Hepa img767.gif i.e. Hepa img768.gif .

However, from the numerical point of view the two forms are not equivalent, and second-order equations are most efficiently treated by special methods.

The general solution of a second-order ODE contains 2n arbitrary constants, which have to be fixed, e.g. by fixing initial values y(x0)=y0 and Hepa img769.gif at one given x = x0, or by fixing boundary values y(x0)=y0, y(x1)=y1 at two points Hepa img770.gif . For numerical methods for initial and boundary value problems see Hall76, Press95.

An example is the Lorentz equation of motion describing the movement of charged particles in a magnetic field:

Hepa img771.gif

with s the path length along the track, and Hepa img772.gif the momentum, along the direction of the track. B is the magnetic field vector. In the Monte Carlo simulation of tracks, one has to solve an initial value problem. In track reconstruction, one has to determine the initial values y(x0) and Hepa img773.gif from a number of measured values of Hepa img774.gif and x along the track, and this is more like a boundary value problem (we have assumed here that the field B is along z). Required here is an integration method for second-order equations. The bending of tracks often being small, one can get good precision using a high (e.g. fourth) order method with quite large steps.

A typical spectrometer magnet has a very sharp-edged field. For the equation of motion Hepa img765.gif this means that Hepa img759.gif resembles a step function. Certain methods (like n-step methods with n>2 and large steps) do not handle such functions very well. On a smaller scale, Hepa img759.gif may have artificial discontinuities due to a discontinuous representation of the magnetic field, or Hepa img775.gif may be discontinuous. Such discontinuities typically invalidate error estimates, and may cause trouble for methods based on extrapolation to the limit of zero step length.

Runge-Kutta methods are simple and efficient, and are much used for this problem. An interesting alternative is offered by the predictor-corrector methods.