## Solving differential equations numerically

For the numerical calculation (or rather, approximation) of the solution of a differential equation, one uses algorithms that are related to the Euler algorithm. Briefly, the Euler algorithm boils down to the following. Even though one does not know the solution of the differential equation explicitly, one does know the derivative of the solution in every point (this is given by the right hand side of the differential equation), or, in other words, one knows what the direction field of the differential equation is. This direction field is now used to approximate the solution. From a given initial value, the next point of the solution is calculated by moving from the initial value in the direction of the direction field. From the point obtained in this way, one does the same again, etc. So the quality of the approximation of the solution depends on the initial value, the direction field, and the step size.

In MATLAB, differential equations can be solved numerically with the commands ode45, ode23 or ode15s. The underlying algorithms of ode45 and ode23 make use of Runge-Kutta-Fehlberg integration with variable step size, i.e., the algorithm increases the step size when the solution varies less. ode23 uses the second and third order formulas, while ode45 uses the fourth and fifth order formulas. The routine ode15s uses another integration routine. This routine is specialized in problems, that the previous routines have problems to deal with: so-called stiff differential equations.

The mentioned routines can be applied to sets of first order differential equations of the form:

 (5.5)

Here is the state vector and the vector valued function that gives the time-derivative of the state vector as a function of and .

Throughout this Section, extensive use is made of function-files to facilitate numerical solving of differential equation. For an explanation of these, see Section 1.6, under "User defined functions".

The form in which ode23 is used, is:

>> [t,y] = ode23('filename',[t0,t1],y0)

Here `filename' is the name of the .m file in which the differential equation is defined, the time-interval over which the differential equation is solved, and the vector with initial values. The algorithm in the program ode23 first determines in a point the derivative of the function defined in the function file 'filename'. With this, the next point is determined, etc. The algorithm itself determines with which time steps it goes through the interval . Since ode45 works with higher order formulas, less integration steps are needed with this command. As a consequence of this, ode23 in general gives less smooth figures than ode45. The output of ode gives a vector with the time instances at which the solution has been calculated, and a row vector with associated solution values.

Choice of solver
For most differential equations, both ode23 and ode45 are suited. The difference between ode23 and ode45 is rather subtle: with the same accuracy, ode23 will need more time steps, though each time step is calculated faster.

However, for a certain class of differential equations, the stiff differential equations, the previously mentioned solvers will not be able to find an accurate solution, or may need excessive computation times for taking very small time steps. In that case, the problem might be of a class called "Stiff differential equations". For these differential equations, the ode15s is a better choice.

Esteur 2010-03-22