Shooting Methods

Consider a two-point boundary problem that is to be solved from an initial independent variable value x = x to a far-field location x = x _e . Not all of the dependent variables will be known at x = x and some boundary conditions will need to be maintained at x = x _e . One way to solve this type of problem is by using a technique known as shooting. An initial guess is made for the free variables at x = x . The ODE is then integrated out to x = x _e . If the far-field boundary conditions are not met, the x = x free variables are updated. The iteration continues until the x = x _e boundary conditions are met to within a specified tolerance.

The trick with shooting methods is to develop a rapidly converging algorithm to obtain updates to the x = x free variables. The method we will use is called multi-dimensional, globally convergent Newton-Raphson. Consider two arrays, V[] that contains the values of the free variables at the x = x boundary and E[] that contains the difference between the computed value of the boundary condition variables at the x = x _e boundary and their true boundary condition values. The size of the V[] and E[] arrays must be the same. To find an array d V[] that will zero the elements of the E[] array requires the solution of the system of equations shown in Eq. (20.21).

Equation 20.21

The updated values of the free variables at the x = x boundary can then be obtained, Eq. (20.22).

Equation 20.22

Unfortunately, there is no general analytic expression for the dE / dV matrix. The matrix elements can be estimated in a finite-difference manner, as in Eq. (20.23).

Equation 20.23

The process to calculate the D E / D V matrix takes a number of steps. With an initial guess for the free variables at the x = x boundary, integrate the ODE to get an initial value for the error vector at x = x _e . You then increment the first free variable at x = x by a small amount and reintegrate the ODE to see how the E[] array values change. Subtracting the updated E[] array from the original gives you the first column of the D E / D V matrix. You then continue the process by incrementing the other free variables at x = x and reintegrating the ODE until the entire D E / D V matrix is filled. You then invert the matrix and obtain the updates to the V[] array.

It may take more than one update to the V[] array until the far-field boundary conditions are satisfied. This sounds like a lot of work and it is when compared to the solution process for initial value problems. However, with a reasonably good initial guess for the free variables at x = x , the solution should converge within three to four iterations.

One final conceptual note about shooting is that sometimes if the initial guess for V[] is not very good the computed updates d V may overshoot physically allowable values. Depending on the ODE, you might get things like a square root of a negative number when you evaluate the ODE right-hand side. One way to enhance the stability of the solution process is to scale the d V updates by a number between 0 and 1. This procedure is called under-relaxation .

We will write an ODEshooter() method that will implement the ODE shooting technique we have just described. This method, like the others that preceded it in this chapter, will be placed inside the ODESolver class and will be a public , static method. The ODEshooter() method makes use of the matrix inversion method EqnSolver.invertMatrix() from Chapter 19 so an appropriate import declaration is placed at the top of the ODESolver class code listing.

After initializing the dependent variables at x = x , the ODEshooter() method solves the ODE by calling the embeddedRK5() method. A first calculation of the error at the x = x _e boundary is performed by having the ODE object call its getError() method. Any ODE subclass that represents a two-point boundary problem will override the getError() method from the ODE class to compute error in the proper manner for that ODE.

The ODEshooter() method then enters a while() loop that updates the free variables at x = x until the error in the far-field boundary conditions is below a specified tolerance. The V[] array updates are under- relaxed by a factor of 0.5. When convergence is achieved, the method exits and returns the number of dependent variable steps used to integrate the ODE. The ODEshooter() method source code is shown next .

 public static int ODEshooter(ODE ode, double V[],           double range, double dx, double tolerance) {   //  Define some convenience variables to make   //  the code more readable.   int numEqns = ode.getNumEqns();   int numVar = ode.getNumFreeVariables();   double x[] = ode.getX();   double y[][] = ode.getY();   //  define some local variables. The E[]  array   //  holds the error at the end of the range of   //  integration.   double E[] = new double[numVar];   double dxInit = dx;   double maxE, dVtotal;   double deltaV = 0.0001;   double underRelax = 0.5;   int i, j, numSteps;   double dV[][] = new double[numVar][numVar];   double dEdV[][] = new double[numVar][numVar];   double Etmp[] = new double[numVar];   //  load initial conditions   ode.setInitialConditions(V);   //  Solve the ODE over the desired range and compute   //  the initial error at the end of the range   numSteps =       ODESolver.embeddedRK5(ode, range, dx, tolerance);   ode.getError(E, y[numSteps-1]);   //  If the E[]  array doesn't meet the desired   //  tolerance try again with new initial conditions.   maxE = 0.0;   for(i=0; i<numVar; ++i) {     maxE = Math.max(maxE,Math.abs(E[i]));   }   while (maxE > tolerance) {     //  Fill the dV array.  Each row of the array     //  is the original V  array with one of its     //  elements perturbed.     for(i=0; i<numVar; ++i) {       for(j=0; j<numVar; ++j) {         dV[i][j] = V[j];       }       dV[i][i] += deltaV;     }     //  Fill the dEdV matrix by determining how the E[]     //  elements change when one of the V[] elements is     //  incremented.     for(j=0; j<numVar; ++j) {       //  Set initial conditions for a given row       ode.setInitialConditions(dV[j]);       dx = dxInit;       //  Solve ODE again with V+dVj       numSteps =        ODESolver.embeddedRK5(ode, range, dx, tolerance);       //  Recompute error for V+dVj.       ode.getError(Etmp, y[numSteps-1]);       //  Compute dEdV       for(i=0; i<numVar; ++i ) {         dEdV[i][j] = ( Etmp[i] - E[i] )/deltaV;       }     }     //  Invert dEdV matrix     EqnSolver.invertMatrix(dEdV);     //  Update V[] matrix. The updates to V[] are     //  under-relaxed to enhance stability     for(i=0; i<numVar; ++i) {       dVtotal = 0.0;       for(j=0; j<numVar; ++j) {         dVtotal += -dEdV[i][j]*E[j];       }       V[i] += underRelax*dVtotal;     }     //  update initial conditions of y[][] using V[]     ode.setInitialConditions(V);     //  Integrate ODE with new initial conditions     numSteps =        ODESolver.embeddedRK5(ode, range, dx, tolerance);     //  Compute new E[]  array and determine maximum     //  error     ode.getError(E, y[numSteps-1]);     maxE = 0.0;     for(i=0; i<numVar; ++i) {       maxE = Math.max(maxE,Math.abs(E[i]));     }   } //  end of while loop   return numSteps; }