Solving Higher Order, Ordinary Differential Equations

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 how to solve higher order ordinary differential equations. The ordinary differential equations solved in Euler's and 2nd and 4th order Runge-Kutta method files were all of the first order. These methods can be used to solve higher order ordinary differential equations by converting them into a coupled system of first order equations. The following sheet demonstrates how this is done.

Section 1: Input

This differential equation used here to illustrate this method is of the form 
The user may specify values for A1, A2, A3, which may be functions of x and y, and A4 which may only be a function of x.

A1[x_, y_] := 2;
A2[x_, y_] := 4;
A3[x_, y_] := 6;
A4[x_] := 7;

The boundary condition must be specified.

x0 = 0.0;
y0 = 10.0;

The derivateve dy/dx at x = x0

dydx0 = 10.0;

Value of x at which y is desired

xf = 1.0;

Section 2: Exact Solution

soln = 
      DSolve[{A1[x, y[x]]*y''[x] + A2[x, y[x]]*y'[x] + 
              A3[x, y[x]]*y[x] + A4[x] == 0, y'[x0] == dydx0, 
          y[x0] == y0}, y, x]

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

Section 3: Approximate Solution

We are going to use Runge-Kutta's 4th order method of solving ODE's, although any method could be used. Our first task is to make the single second order ODE, a system of two first order ODE's. Thus, the derivative dy/dx is set equal to another variable.

dy/dx = z: This will be function #1, f1 and used for determining k values for the y approximation.

f1[x_, y_, z_] := z;

dz/dx = -(A2*z + A3*y + A4) / A1: This will be function #2, f2 and used for determining k values for the z approximation.

f2[x_, y_, 
          z_] := -1.0/
            A1[x, y]*(A2[x, y]*z + A3[x, y]*y + A4[x]);

The step size will  be used throughout this sheet

h = (xf - x0)/5.0

The following shows five steps of the 4th order Runge-Kutta method for solving  a system of ODEs.

Step 1

Using the values of the initial conditions, we find the value of the solution one step forward.

X[0] = x0
Y[0] = y0
Z[0] = dydx0

X[1] = X[0] + h

Here the k values are calculated. The ky values use function #1, f1 for estimates and the kz values use function #2, f2. Note how each k values is based in part on the preceding k value.

k1y = f1[X[0], Y[0], Z[0]]

k1z = f2[X[0], Y[0], Z[0]]

k2y = f1[X[0] + 0.5*h, Y[0] + 0.5*k1y*h, Z[0] + 0.5*k1z*h]

k2z = f2[X[0] + 0.5*h, Y[0] + 0.5*k1y*h, Z[0] + 0.5*k1z*h]

k3y = f1[X[0] + 0.5*h, Y[0] + 0.5*k2y*h, Z[0] + 0.5*k2z*h]

k3z = f2[X[0] + 0.5*h, Y[0] + 0.5*k2y*h, Z[0] + 0.5*k2z*h]

k4y = f1[X[0] + h, Y[0] + k3y*h, Z[0] + k3z*h]

k4z = f2[X[0] + h, Y[0] + k3y*h, Z[0] + k3z*h]

With all the k values calculated we can solve for the solution at the next step and the derivative (used in next step calculation.)

Y[1] = Y[0] + 1/6*(k1y + 2*k2y + 2*k3y + k4y)*h

Z[1] = Z[0] + 1/6*(k1z + 2*k2z + 2*k3z + k4z)*h

Step 2

X[2] = X[1] + h

k1y = f1[X[1], Y[1], Z[1]]

k1z = f2[X[1], Y[1], Z[1]]

k2y = f1[X[1] + 0.5*h, Y[1] + 0.5*k1y*h, Z[1] + 0.5*k1z*h]

k2z = f2[X[1] + 0.5*h, Y[1] + 0.5*k1y*h, Z[1] + 0.5*k1z*h]

k3y = f1[X[1] + 0.5*h, Y[1] + 0.5*k2y*h, Z[1] + 0.5*k2z*h]

k3z = f2[X[1] + 0.5*h, Y[1] + 0.5*k2y*h, Z[1] + 0.5*k2z*h]

k4y = f1[X[1] + h, Y[1] + k3y*h, Z[1] + k3z*h]

k4z = f2[X[1] + h, Y[1] + k3y*h, Z[1] + k3z*h]

Y[2] = Y[1] + 1/6*(k1y + 2*k2y + 2*k3y + k4y)*h

Z[2] = Z[1] + 1/6*(k1z + 2*k2z + 2*k3z + k4z)*h

Step 3

X[3] = X[2] + h

k1y = f1[X[2], Y[2], Z[2]]

k1z = f2[X[2], Y[2], Z[2]]

k2y = f1[X[2] + 0.5*h, Y[2] + 0.5*k1y*h, Z[2] + 0.5*k1z*h]

k2z = f2[X[2] + 0.5*h, Y[2] + 0.5*k1y*h, Z[2] + 0.5*k1z*h]

k3y = f1[X[2] + 0.5*h, Y[2] + 0.5*k2y*h, Z[2] + 0.5*k2z*h]

k3z = f2[X[2] + 0.5*h, Y[2] + 0.5*k2y*h, Z[2] + 0.5*k2z*h]

k4y = f1[X[2] + h, Y[2] + k3y*h, Z[2] + k3z*h]

k4z = f2[X[2] + h, Y[2] + k3y*h, Z[2] + k3z*h]

Y[3] = Y[2] + 1/6*(k1y + 2*k2y + 2*k3y + k4y)*h

Z[3] = Z[2] + 1/6*(k1z + 2*k2z + 2*k3z + k4z)*h

Step 4

X[4] = X[3] + h

k1y = f1[X[3], Y[3], Z[3]]

k1z = f2[X[3], Y[3], Z[3]]

k2y = f1[X[3] + 0.5*h, Y[3] + 0.5*k1y*h, Z[3] + 0.5*k1z*h]

k2z = f2[X[3] + 0.5*h, Y[3] + 0.5*k1y*h, Z[3] + 0.5*k1z*h]

k3y = f1[X[3] + 0.5*h, Y[3] + 0.5*k2y*h, Z[3] + 0.5*k2z*h]

k3z = f2[X[3] + 0.5*h, Y[3] + 0.5*k2y*h, Z[3] + 0.5*k2z*h]

k4y = f1[X[3] + h, Y[3] + k3y*h, Z[3] + k3z*h]

k4z = f2[X[3] + h, Y[3] + k3y*h, Z[3] + k3z*h]

Y[4] = Y[3] + 1/6*(k1y + 2*k2y + 2*k3y + k4y)*h

Z[4] = Z[3] + 1/6*(k1z + 2*k2z + 2*k3z + k4z)*h

Step 5

X[5] = X[4] + h

k1y = f1[X[4], Y[4], Z[4]]

k1z = f2[X[4], Y[4], Z[4]]

k2y = f1[X[4] + 0.5*h, Y[4] + 0.5*k1y*h, Z[4] + 0.5*k1z*h]

k2z = f2[X[4] + 0.5*h, Y[4] + 0.5*k1y*h, Z[4] + 0.5*k1z*h]

k3y = f1[X[4] + 0.5*h, Y[4] + 0.5*k2y*h, Z[4] + 0.5*k2z*h]

k3z = f2[X[4] + 0.5*h, Y[4] + 0.5*k2y*h, Z[4] + 0.5*k2z*h]

k4y = f1[X[4] + h, Y[4] + k3y*h, Z[4] + k3z*h]

k4z = f2[X[4] + h, Y[4] + k3y*h, Z[4] + k3z*h]

Y[5] = Y[4] + 1/6*(k1y + 2*k2y + 2*k3y + k4y)*h

Z[5] = Z[4] + 1/6*(k1z + 2*k2z + 2*k3z + k4z)*h

Results

Approximate solution at x = xf

Y[5]

Exact solution at x = xf

yf = Evaluate[y[xf]  /.  soln][[1]]

Absolute Relative True Error Percentage

Abs[(Y[5] - yf)/yf]*100

Graphical depiction of the solution

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

Show[plot1, plot2]

References

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

Conclusions

This sheet has showed us how we can use the methods learned for solving first order ordinary differential equaitons and apply them to solve higher order problems. In this example we just did a second order ordinary differential equation, but the same method could be utilitized to solve higher orders as well.