Shooting Method of Solving 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 the shooting method of solving ordinary differential equations. The shooting method of solving ODEs is used when there are multiple initial conditions to be satisfied at different independent variable locations. For example, as in many real life problems such as a simply supported beam, the boundary conditions are y(0) = 0 and y(L) = 0 where y represents the displacement and L is the beam length. The methods for solving ODEs begin at one point and using approximations of the derivative, find the solution from one end to another. These methods will not take into account the bondary condition on the other end. The shooting method uses a method of solving ODEs (4th order Runge-Kutta in this example) in a way that satisfies both boundary conditions.

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_] := 1;
A2[x_, y_] := 1;
A3[x_, y_] := 1;
A4[x_] := 1;

The boundary condition must be specified.

x0 = 0.0;
y0 = 10.0;

xf = 1.0;
yf = 0.0;

An initial value for the derivative dy/dx  is assumed and will be adjusted to satisfy the second boundary condition.

dydx1 = -15.0;

An second value for the derivative dy/dx is assumed here as well to be used in second iteration.

dydx2 = -10.0;

Section 2: Setup of Approximate Solution

To use Runge-Kutta's 4th order method of solving ODEs, the equations must be of order one. Thus, the derivative dy/dx is set equal to another variable and a system of first order differential equations is made.

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.

Section 3: First Iteration

On the first iteration we will show each calculation of each step. This is so the method can be carefully analyzed and understood.

Step 1

First the initial condition at x = x0 is used to begin the 4th order Runge-Kutta method. The first assumed value of dy/dx at x = x0 is used as well.

X1[0] = x0
Y1[0] = y0
Z1[0] = dydx1

X1[1] = X1[0] + h

The k values are now calculated for y and z.

k1y = f1[X1[0], Y1[0], Z1[0]]

k1z = f2[X1[0], Y1[0], Z1[0]]

k2y = 
      f1[X1[0] + 0.5*h, Y1[0] + 0.5*k1y*h, Z1[0] + 0.5*k1z*h]

k2z = 
      f2[X1[0] + 0.5*h, Y1[0] + 0.5*k1y*h, Z1[0] + 0.5*k1z*h]

k3y = 
      f1[X1[0] + 0.5*h, Y1[0] + 0.5*k2y*h, Z1[0] + 0.5*k2z*h]

k3z = 
      f2[X1[0] + 0.5*h, Y1[0] + 0.5*k2y*h, Z1[0] + 0.5*k2z*h]

k4y = f1[X1[0] + h, Y1[0] + k3y*h, Z1[0] + k3z*h]

k4z = f2[X1[0] + h, Y1[0] + k3y*h, Z1[0] + k3z*h]

Once the k values have been found, y and z can be determined. z will be used in the next iteration.

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

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

Step 2

X1[2] = X1[1] + h

k1y = f1[X1[1], Y1[1], Z1[1]]

k1z = f2[X1[1], Y1[1], Z1[1]]

k2y = 
      f1[X1[1] + 0.5*h, Y1[1] + 0.5*k1y*h, Z1[1] + 0.5*k1z*h]

k2z = 
      f2[X1[1] + 0.5*h, Y1[1] + 0.5*k1y*h, Z1[1] + 0.5*k1z*h]

k3y = 
      f1[X1[1] + 0.5*h, Y1[1] + 0.5*k2y*h, Z1[1] + 0.5*k2z*h]

k3z = 
      f2[X1[1] + 0.5*h, Y1[1] + 0.5*k2y*h, Z1[1] + 0.5*k2z*h]

k4y = f1[X1[1] + h, Y1[1] + k3y*h, Z1[1] + k3z*h]

k4z = f2[X1[1] + h, Y1[1] + k3y*h, Z1[1] + k3z*h]

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

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

Step 3

X1[3] = X1[2] + h

k1y = f1[X1[2], Y1[2], Z1[2]]

k1z = f2[X1[2], Y1[2], Z1[2]]

k2y = 
      f1[X1[2] + 0.5*h, Y1[2] + 0.5*k1y*h, Z1[2] + 0.5*k1z*h]

k2z = 
      f2[X1[2] + 0.5*h, Y1[2] + 0.5*k1y*h, Z1[2] + 0.5*k1z*h]

k3y = 
      f1[X1[2] + 0.5*h, Y1[2] + 0.5*k2y*h, Z1[2] + 0.5*k2z*h]

k3z = 
      f2[X1[2] + 0.5*h, Y1[2] + 0.5*k2y*h, Z1[2] + 0.5*k2z*h]

k4y = f1[X1[2] + h, Y1[2] + k3y*h, Z1[2] + k3z*h]

k4z = f2[X1[2] + h, Y1[2] + k3y*h, Z1[2] + k3z*h]

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

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

Step 4

X1[4] = X1[3] + h

k1y = f1[X1[3], Y1[3], Z1[3]]

k1z = f2[X1[3], Y1[3], Z1[3]]

k2y = 
      f1[X1[3] + 0.5*h, Y1[3] + 0.5*k1y*h, Z1[3] + 0.5*k1z*h]

k2z = 
      f2[X1[3] + 0.5*h, Y1[3] + 0.5*k1y*h, Z1[3] + 0.5*k1z*h]

k3y = 
      f1[X1[3] + 0.5*h, Y1[3] + 0.5*k2y*h, Z1[3] + 0.5*k2z*h]

k3z = 
      f2[X1[3] + 0.5*h, Y1[3] + 0.5*k2y*h, Z1[3] + 0.5*k2z*h]

k4y = f1[X1[3] + h, Y1[3] + k3y*h, Z1[3] + k3z*h]

k4z = f2[X1[3] + h, Y1[3] + k3y*h, Z1[3] + k3z*h]

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

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

Step 5

X1[5] = X1[4] + h

k1y = f1[X1[4], Y1[4], Z1[4]]

k1z = f2[X1[4], Y1[4], Z1[4]]

k2y = 
      f1[X1[4] + 0.5*h, Y1[4] + 0.5*k1y*h, Z1[4] + 0.5*k1z*h]

k2z = 
      f2[X1[4] + 0.5*h, Y1[4] + 0.5*k1y*h, Z1[4] + 0.5*k1z*h]

k3y = 
      f1[X1[4] + 0.5*h, Y1[4] + 0.5*k2y*h, Z1[4] + 0.5*k2z*h]

k3z = 
      f2[X1[4] + 0.5*h, Y1[4] + 0.5*k2y*h, Z1[4] + 0.5*k2z*h]

k4y = f1[X1[4] + h, Y1[4] + k3y*h, Z1[4] + k3z*h]

k4z = f2[X1[4] + h, Y1[4] + k3y*h, Z1[4] + k3z*h]

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

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

Results

Now that the five steps have been calculated we plot the results and examine how close the value at y1[5] has come to the boundary condition at that value.

plot1 = 
        ListPlot[{{x0, y0}, {xf, 
              yf}}, 
PlotStyle -> (PointSize[
                0.03])];

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

Show[plot1, plot2];

Section 4: Second Iteration

Now the calculation in the first iteration is repeated but for a different assumed dy/dx at x = x0. This calculation is done the same way as before, but here we have enclosed it in a loop to save space.

X2[0] = x0
Y2[0] = y0
Z2[0] = dydx2

For[i = 0, 
      i < 5, i++, 
X2[i + 1] = 
        X2[i] + h; 
k1y = 
        f1[X2[i], Y2[i], Z2[i]]; 
k1z = 
        f2[X2[i], Y2[i], Z2[i]]; 
k2y = 
        f1[X2[i] + 0.5*h, Y2[i] + 0.5*k1y*h, 
          Z2[i] + 0.5*k1z*h]; 
k2z = 
        f2[X2[i] + 0.5*h, Y2[i] + 0.5*k1y*h, 
          Z2[i] + 0.5*k1z*h]; 
k3y = 
        f1[X2[i] + 0.5*h, Y2[i] + 0.5*k2y*h, Z2[i] + 0.5*k2z*h]; 
      k3z = f2[X2[i] + 0.5*h, Y2[i] + 0.5*k2y*h, 
          Z2[i] + 0.5*k2z*h]; 
k4y = 
        f1[X2[i] + h, Y2[i] + k3y*h, Z2[i] + k3z*h]; 
k4z = 
        f2[X2[i] + h, Y2[i] + k3y*h, Z2[i] + k3z*h]; 
Y2[
          i + 1] = 
        Y2[i] + 1/6*(k1y + 2*k2y + 2*k3y + k4y)*h; 
Z2[
          i + 1] = 
        Z2[i] + 1/6*(k1z + 2*k2z + 2*k3z + k4z)*
            h
]

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

Show[plot1, plot2, plot3];

Section 5: Third Iteration

By examining the figure immediately above, you will note that at x = x0, the solution (for both iterations) is close to the boundary condition. At x = xf  however, it is likely that the solution does not (unless you are a good guesser). We can use the first two iterations to predict what a better value for dy/dx should be. We will assume that the relationship is linear and use interpolation to determine the final estimate for dy/dx. This will yield exact results if the equation is linear while for nonlinear differential equations it will be a reasonable approximation.

The third assumed dy/dx.

dydx3 = (dydx1 - dydx2)/(Y1[5] - Y2[5])*(yf - Y2[5]) + 
        dydx2

X3[0] = x0
Y3[0] = y0
Z3[0] = dydx3

As in the second iteration we are doing this calculation in a loop.

For[i = 0, 
      i < 5, i++, 
X3[i + 1] = 
        X3[i] + h; 
k1y = 
        f1[X3[i], Y3[i], Z3[i]]; 
k1z = 
        f2[X3[i], Y3[i], Z3[i]]; 
k2y = 
        f1[X3[i] + 0.5*h, Y3[i] + 0.5*k1y*h, 
          Z3[i] + 0.5*k1z*h]; 
k2z = 
        f2[X3[i] + 0.5*h, Y3[i] + 0.5*k1y*h, 
          Z3[i] + 0.5*k1z*h]; 
k3y = 
        f1[X3[i] + 0.5*h, Y3[i] + 0.5*k2y*h, Z3[i] + 0.5*k2z*h]; 
      k3z = f2[X3[i] + 0.5*h, Y3[i] + 0.5*k2y*h, 
          Z3[i] + 0.5*k2z*h]; 
k4y = 
        f1[X3[i] + h, Y3[i] + k3y*h, Z3[i] + k3z*h]; 
k4z = 
        f2[X3[i] + h, Y3[i] + k3y*h, Z3[i] + k3z*h]; 
Y3[
          i + 1] = 
        Y3[i] + 1/6*(k1y + 2*k2y + 2*k3y + k4y)*h; 
Z3[
          i + 1] = 
        Z3[i] + 1/6*(k1z + 2*k2z + 2*k3z + k4z)*
            h
]

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

Show[plot1, plot2, plot3, plot4];