Runge-Kutta 4th Order Method - Convergence

Nathan Collier, Autar Kaw
University of South Florida
United States of America
kaw@eng.usf.edu

Introduction

This worksheet demonstrates the use of Mathematica to illustrate the convergence of the Runge-Kutta 4th order method of solving ordinary differential equations. This method of solving ordinary differential equations uses the derivative and value of the function at the boundary condition to project the location and value of the next point on the function.

Section 1: Input

The following simulation illustrates the convergence of the Runge-Kutta 4th order method of solving ordinary differential equations (ODEs). This section is the only section where the user interacts with the program.  The user enters ordinary differential equation of the form  dy/dx = f(x, y), the boundary condition, and the value of x at which the solution is desired. By entering this data, the program will calculate the exact (Mathematica numerical value if it is not exact) value of the solution, followed by the results using the method with 1, 2, 4, 8 ... n steps. The program will also display the true error, the absolute relative percentage true error, the approximate error, the absolute relative approximate percentage error, and the least number of significant digits correct all as a function of number of segments.

Ordinary differential equation of the form dy/dx = f(x,y)

f[x_, y_] := y*x^2 - 1.2*y;

Boundary condition for x

x0 = 0.0;

Boundary condition for y

y0 = 1.0;

Value of x at which y is desired

xf = 2.0;

Maximum number of steps, n. Valid entries are 1, 2, 4, 8 ... 128.

n = 128;

Section 2: Procedure

The following procedure estimates the solution of ordinary differential equations at a point xf.
n = number of steps
x0 = boundary condition for x
y0 = boundary condition for  y
xf = value of x at which y is desired
f = differential equation in the form dy/dx = f(x,y)

RK4th[n_, x0_, y0_, xf_, f_] := 
      Module[{Y, h, i, k1, k2, k3, 
          k4}, 
h = (xf - x0)/n; 
Y = 
          y0; 
For[i = 0, 
          i < n, i++, 
k1 = 
            f[x0 + i*h, Y]; 
k2 = 
            f[x0 + i*h + 0.5*h, Y + 0.5*k1*h]; 
k3 = 
            f[x0 + i*h + 0.5*h, Y + 0.5*k2*h]; 
k4 = 
            f[x0 + i*h + h, Y + k3*h]; 
Y = 
            Y + 1/6*(k1 + 2*k2 + 2*k3 + k4)*
                h]; 
Y]

Section 3: Calculation

The exact value of the solution of the ODE (EV):

soln = 
      DSolve[{f[x, y[x]] == y'[x], y[x0] == y0}, y, x]

EV = (y /. First[soln])[xf]

plot1 = 
        Plot[Evaluate[y[x]  /.  soln], {x, x0, xf}, 
          PlotStyle -> {{Thickness[0.012], 
                RGBColor[1, 0, 0]}}];

This loop calculates the following
AV = approximate value of the ODE solution using Euler's Method by calling the module EulerMethod
Et = true error
et = absolute relative true error percentage
Ea = approximate error
ea = absolute relative approximate error percentage
sig = least number of significant digits correct in approximation

nnnn = Log[n]/Log[2];

nth = 7;
Do[
Nn[i] = 
          2^i; 
H[i] = (xf - x0)/
            Nn[i]; 
AV[i] = 
          RK4th[2^i, x0, y0, xf, f]; 
Et[i] = 
          EV - AV[i]; 
et[i] = 
          Abs[(Et[i]/EV)]*100.0; 
If[
          i > 0, 
Ea[i] = 
            AV[i] - AV[i - 1]; 
ea[i] = 
            Abs[Ea[i]/AV[i]]*100.0; 
sig[i] = 
            Floor[(2 - Log[10, ea[i]/0.5])]; 
If[
            sig[i] < 0, 
            sig[i] = 0];
]
, {i, 0, 
          nth}];

Section 5: Graphs

The following code uses the method to calculate intermediate step values for the purpose of displaying the solution over the entire range specified. The number of steps used is the maximum value, n, defined in the Input section.

X[0] = x0;
Y[0] = y0;
h = (xf - x0)/n;
For[i = 1, 
        i <= n, i++, 
X[i] = 
          x0 + i*h; 
k1 = 
          f[X[i - 1], Y[i - 1]]; 
k2 = 
          f[X[i - 1] + 0.5*h, Y[i - 1] + 0.5*k1*h]; 
k3 = 
          f[X[i - 1] + 0.5*h, Y[i - 1] + 0.5*k2*h]; 
k4 = 
          f[X[i - 1] + h, Y[i - 1] + k3*h]; 
Y[i] = 
          Y[i - 1] + 1/6*(k1 + 2*k2 + 2*k3 + k4)*h];

data = Table[{X[i], Y[i]}, {i, 0, n}];
plot2 = 
        ListPlot[
          data, 
PlotJoined -> 
            True, 
PlotStyle -> {Thickness[0.012], 
              RGBColor[0, 1, 0]}, 
DisplayFunction -> 
            Identity];

Show[plot1, plot2,  
      PlotLabel -> "Exact and Approximate Solution"]

data = 
        Table[{H[i], AV[i]}, {i, 0, nth}];
\(\(plot2 = 
        ListPlot[
          data, 
PlotJoined -> 
            True, 
PlotStyle -> {Thickness[0.012], 
              RGBColor[1, 0, 
                0]}, 
PlotLabel -> "Approximate value of the solution of the ODE
at x = xf as a function of step size"];\)\)

data = 
        Table[{Nn[i], AV[i]}, {i, 0, nth}];
plot2 = 
        ListPlot[
          data, 
PlotJoined -> 
            True, 
PlotStyle -> {Thickness[0.012], 
              RGBColor[1, 0, 
                0]}, 
PlotLabel -> "Approximate value of \
the solution of the ODE\nat x = xf as a function of number of \
steps"];

data = 
        Table[{Nn[i], Et[i]}, {i, 0, nth}];
plot2 = 
        ListPlot[
          data, 
AxesOrigin -> {0, 
              0}, 
PlotJoined -> 
            True, 
PlotStyle -> {Thickness[0.012], 
              RGBColor[0, 0, 
                1]}, 
PlotLabel -> "True error as a \
function of number of steps"];

data = 
        Table[{Nn[i], et[i]}, {i, 0, nth}];
plot2 = 
        ListPlot[
          data, 
AxesOrigin -> {0, 
              0}, 
PlotJoined -> 
            True, 
PlotStyle -> {Thickness[0.012], 
              RGBColor[0, 0, 
                1]}, 
PlotLabel -> "Absolute relative \
true error percentage \nas a function of number of steps"];

data = 
        Table[{Nn[i], Ea[i]}, {i, 1, nth}];
plot2 = 
        ListPlot[
          data, 
AxesOrigin -> {0, 
              0}, 
PlotJoined -> 
            True, 
PlotStyle -> {Thickness[0.012], 
              RGBColor[0, 1, 
                0]}, 
PlotLabel -> "Approximate error \n\
as a function of number of steps"];

data = 
        Table[{Nn[i], ea[i]}, {i, 1, nth}];
plot2 = 
        ListPlot[
          data, 
AxesOrigin -> {0, 
              0}, 
PlotJoined -> 
            True, 
PlotStyle -> {Thickness[0.012], 
              RGBColor[0, 1, 
                0]}, 
PlotLabel -> "Absolute relative \
approximate error percentage \nas a function of number of steps"];

data = 
        Table[{Nn[i], sig[i]}, {i, 1, nth}];
plot2 = 
        ListPlot[
          data, 
AxesOrigin -> {0, 
              0}, 
PlotJoined -> 
            True, 
PlotStyle -> {Thickness[0.012], 
              RGBColor[0.80, 0.52, 
                0.25]}, 
PlotLabel -> "Least correct \
significant digits \nas a function of number of steps"];

References

Chapra & Canale, Numerical Methods for Engineering: With Programming Applications. 3rd ed. WCB/McGraw-Hill. 1998.

Conclusions

Mathcad helped to see that the Runge-Kutta 4th order method of solving ODEs converges very quickly.