8th Order Runge-Kutta for Integrating System of ODEs (ODE86)

ODE86 Integrates a system of ordinary differential equations using  a 12-stages, 8th and 6th order Runge-Kutta formulas. Better than ODE45 for tolerances stringent than 1e-6.

Forms of using the function ODE86 :

• [t,y] = ODE86(‘yprime’, t0, tfinal, y0) – integrates the system of ordinary differential equations described by the M-file main.m over the interval T0 to Tfinal and using initial conditions Y0.
• [t, y] = ODE86(f, t0, tfinal, y0, tol) uses tolerance TOL.

INPUT:

• f – String containing name of user-supplied problem description.
• t – Time (scalar).
• y – Solution column-vector.
• yprime – Returned derivative column-vector; yprime(i) = dy(i)/dt.
• t0 – Initial value of t.
• tfinal- Final value of t.
• y0 – Initial value column-vector.
• tol – The desired accuracy. (Default: tol = 1.e-8).

OUTPUT:

• t – Returned integration time points (row-vector).
• y – Returned solution, one solution column-vector per tout-value.
• ir – Returned number of rejected steps

The result can be displayed by: plot(tout, yout). Example: Solve two-body problem using inline function:
The problem :
y1’=y3, y2’=y4, y3’=-y1/(y1^2+y2^2)^1.5, y4’=-y2/(y1^2+y2^2)^1.5
With initial condtitions y1(0)=.5, y2(0)=0, y3(0)=0, y4(0)=3^0.5
Matlab call :
[x,y]=ode86(inline(‘[y(3);y(4);-y(1)/sqrt(y(1)^2+y(2)^2)^3;-y(2)/sqrt(y(1)^2+y(2)^2)^3]’,’x’,’y’), 0, 20, [.5 0 0 sqrt(3)]’, 1e-11);

The code based on the ODE45 by C.B. Moler, 25-3-1987.

The error control method and coefficients are taken from Ch. Tsitouras and S. N. Papakostas, ‘Cheap Error Estimation for Runge-Kutta methods’, SIAM J. Sci.