The Dormand–Prince method for ODEs Version 1

[Namir]

The Dormand–Prince method for ODEs Version 1

The Dormand–Prince method is more accurate than the Runge–Kutta–Fehlberg method and it is used by the MATLAB ode45 solver.

NOTE: Make sure you supply the ODE function with a clear matrix that returns the columns of x and y values.

The parameters of function ODE are:

1. The parameters x0 and y0 represent the initial integration point.
2. The parameter xf is the final x value where integratin stops.
3. The parameter h is the integration step.
4. The parameter hmin is the minimum integration step.
5. The parameter hReduceFactor is the reducion factor for the integration step h.
6. The parameter toler is the error tolerance for the variable y.
7. The parameter resMat is the matrix that returns the columns of sampled x and y values.
8. The parameter skipSampling specifies the number of (x, y) point to exclude from appearing in the resMat matrix. To include all point supply an argument of zero. To skip a sequence of n points supply an argument of n. The function uses an IF statement to ensure that (xf, y) is included in the resMat matrix.

The function returns the matrix resMat populated with (x, y) values.

Example
=======

To integrate fx(x,y) = x+y from (1,2) to xf = 2 in steps of 0.01, minimum step of 0.001, step reduction factor of 2, tolerance of 1e-6, and sampling all values of (x, y):

M1 := ODE(1,2,2,0.01,0.001,2,1e-6,M1,0)

The program creates a matrix with x,y data. The last value of y at x=2 is calculate as 7.87312731391.

Reference
=========

Mathematical Modeling in Chemical Engineering, by Anders Rasmuson, Bengt Andersson, Louise Olsson, and Ronnie Andersson, published by Cambridge University Press 2014. See page 89.

Program Listing
===============

Code:

EXPORT fx(x,y)
BEGIN
  // INTEGRAL IS X^2/2 + X*Y + const
  RETURN x+y;
END;

EXPORT ODE(x0,y0,xf,h,hmin,hReduceFactor,toler,resMat,skipSampling)
BEGIN
  LOCAL i, j, z;
  LOCAL k1, k2, k3, k4, k5, k6, k7;
  LOCAL err, x, y, dy, ok, counter;
  
  counter := skipSampling;
  x := x0;
  y := y0;
  i := 1;
  resMat:=MAKEMAT(0,1,1);  
  resMat(i,1) := x;
  resMat(i,2) := y;
  WHILE x < xf DO
    k1 := fx(x,y);
    k2 := fx(x+h/5, y+h/5*k1);
    k3 := fx(x+3/10*h, y+h*3/40*(k1+3*k2));
    k4 := fx(x+4/5*h, y+h*(44/45*k1-56/15*k2+32/9*k3));
    k5 := fx(x+8/9*h, y+h*(19372/6561*k1-25360/2187*k2+64448/6561*k3-212/729*k4));
    k6 := fx(x+h,y+h*(9017/3168*k1-355/33*k2+46732/5247*k3+49/176*k4-5103/18656*k5));
    z := y + h*(35/384*k1+500/1113*k3+125/192*k4-2187/6784*k5+11/84*k6);
    k7 := fx(x+h, z);
    dy := h*(5179/57600*k1+7571/16695*k3+393/640*k4-92097/339200*k5+187/2100*k6+k7/40);
    err := h*(71/57600*k1-71/16695*k3+71/1920*k4-17253/339200*k5+22/525*k6-k7/40);
    ok := 1;
    IF ABS(err) > toler THEN
      ok := 0;
      IF h / hReduceFactor >= hmin THEN
        h := h / hReduceFactor;
      ELSE
        h := hmin;
        ok := 1;
      END;
    END;
    IF ok == 1 THEN
      x := x + h;
      y := y + dy;
      // sample this point?
      IF counter < 1 THEN
        i := i + 1;
        resMat(i,1) := x;
        resMat(i,2) := y;
        counter := skipSampling;
      ELSE
        counter := counter - 1;
      END;
    END;
  END;
  // make sure the ending point is included in the
  // result matrix
  IF resMat(i,1) < xf THEN
    i := i + 1;
    resMat(i,1) := x;
    resMat(i,2) := y;
  END;
  RETURN resMat;
END;