The Dormand–Prince method for ODEs Version 2

[Namir]

The Dormand–Prince method for ODEs Version 2

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.

NOTES2: This version will alter the integation steps between a user-specified range of steps, given a range of min/max tolernce error values for y.

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 hMinMax is the array of minimum and maximum integration steps.
5. The parameter hReduceFactor is the reducion factor for the integration step h.
6. The parameter tolerMinMax is the array of minimu and maximum error tolerances 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 fx2(x,y) = x+y from (1,2) to xf = 2 in steps of 0.01, minimum step of 0.001, maximum step of 0.05, step reduction factor of 2, minimum tolerance of 1e-6, maximum tolernce of 1e-4, and sampling all values of (x, y):

M1 := ODE2(1,2,2,0.01,[0.001,0.05],2,[1e-6,1e-4],M1,0)

The program creates a matrix with the x, y values. The last value of y at x=2 is calculate as 7.79079503593.

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 fx2(x,y)
BEGIN
  // INTEGRAL IS X^2/2 + X*Y + const
  RETURN x+y;
END;

EXPORT ODE2(x0,y0,xf,h,hMinMax,hReduceFactor,tolerMinMax,resMat,skipSampling)
BEGIN
  LOCAL i, j, z;
  LOCAL hmin, hmax, tolerMin, tolerMax;
  LOCAL k1, k2, k3, k4, k5, k6, k7;
  LOCAL err, x, y, dy, ok, counter;
  
  hmin := hMinMax(1);
  hmax := hMinMax(2);
  tolerMin := tolerMinMax(1);
  tolerMax := tolerMinMax(2);
  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 := fx2(x,y);
    k2 := fx2(x+h/5, y+h/5*k1);
    k3 := fx2(x+3/10*h, y+h*3/40*(k1+3*k2));
    k4 := fx2(x+4/5*h, y+h*(44/45*k1-56/15*k2+32/9*k3));
    k5 := fx2(x+8/9*h, y+h*(19372/6561*k1-25360/2187*k2+64448/6561*k3-212/729*k4));
    k6 := fx2(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 := fx2(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) > tolerMin THEN
      ok := 0;
      IF h / hReduceFactor >= hmin THEN
        h := h / hReduceFactor;
      ELSE
        h := hmin;
        ok := 1;
      END;
    END;
    IF ABS(err) < tolerMax THEN
      ok := 0;
      IF h * hReduceFactor <= hmax THEN
        h := h * hReduceFactor;
      ELSE
        h := hmax;
        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;