Embedded Runge-Kutta Solvers


Embedded Runge -Kutta Solvers

One problem with the fourth-order Runge-Kutta method we developed in the previous section is that it uses a constant independent variable increment, D x , over the entire integration range. In certain situations, regions of high gradients may require smaller step sizes while regions of lesser gradients can maintain solution accuracy with larger step sizes. Using the high-gradient step size over the entire integration range results in an inefficient algorithm. Using a low- gradient step size over the entire range may cause the solution to be inaccurate in certain regions of the domain.

A solution to this problem is to use what is known as an embedded Runge-Kutta technique with adaptive step size control. There are certain types of Runge-Kutta schemes that have embedded within them a lower-order Runge-Kutta scheme. These are called embedded Runge-Kutta algorithms. Why is this significant? Because the difference in solution between the higher- and lower-order schemes can provide an estimate of the truncation error of the solution.

The ability to estimate truncation error is what makes adaptive step size control possible. You can automatically adjust the local step size either up or down so the computed truncation error is within a certain range. There is no longer a need to estimate an appropriate step size for the entire domain. The step size can be quite large in smoothly varying regions of the ODE solution and can be reduced to smaller values in regions of strong gradients.

The embedded Runge-Kutta solver we will use implements the fifth-order Runge-Kutta algorithm in Eq. (20.17).

Equation 20.17

graphics/20equ17.gif


In Eq. (20.17), the a, b, c , and d values are constant coefficients. Using different coefficients, you can also write a fourth-order Runge-Kutta solver, shown in Eq. (20.18), using the same computed D y values.

Equation 20.18

graphics/20equ18.gif


The truncation error at any step in the integration process can be estimated by Eq. (20.19).

Equation 20.19

graphics/20equ19.gif


The a, b, c , and d coefficients we will use were derived by Cash and Karp 1 and are shown in Table 20.1.

The only thing that remains to be done in the development of our embedded Runge-Kutta solver is to come up with a scheme to compute the value of D x for a given step that will keep the truncation error below a specified maximum value. There are several ways to specify the maximum allowable error. You can use a constant value or evaluate the maximum allowable error as a function of the derivatives of the dependent variables . We're going to keep things simple in this example. Since the Runge-Kutta scheme is fifth-order accurate, the error should scale with D x 5 . The optimum step size can be determined from Eq. (20.20).

Table 20.1. Cash-Karp Coefficients

I

A I

B I1

B I2

B I3

B I4

B I5

C I

D I

1

graphics/20inl22.gif

graphics/20inl26.gif

2

graphics/20inl03.gif

graphics/20inl07.gif

3

graphics/20inl04.gif

graphics/20inl08.gif

graphics/20inl12.gif

graphics/20inl23.gif

graphics/20inl27.gif

4

graphics/20inl05.gif

graphics/20inl09.gif

graphics/20inl13.gif

graphics/20inl16.gif

graphics/20inl24.gif

graphics/20inl28.gif

5

1

graphics/20inl10.gif

graphics/20inl14.gif

graphics/20inl17.gif

graphics/20inl19.gif

graphics/20inl29.gif

6

graphics/20inl06.gif

graphics/20inl11.gif

graphics/20inl15.gif

graphics/20inl18.gif

graphics/20inl20.gif

graphics/20inl21.gif

graphics/20inl25.gif

graphics/20inl30.gif

Equation 20.20

graphics/20equ20.gif


The E max parameter is the maximum error or tolerance that you are willing to accept in your calculation. For a situation with more than one dependent variable, Ecurrent would be the maximum current error among the dependent variables.

We are now ready to write a method that implements an embedded Runge-Kutta solver. The method is called embeddedRK5() and is placed inside the ODESolver class previously described in this chapter. In many respects it is similar to the rungeKutta4() method. The embeddedRK5() method takes four arguments, an ODE object and three variables of type double defining the range, initial D x , and error tolerance for the computation.

One principal difference between the embeddedRK5() and rungeKutta4() methods is that the embeddedRK5() method will estimate the maximum truncation error at each step in the integration. If the maximum error is greater than the tolerance, then D x is decreased and the current step is integrated again. If the maximum error is less than the tolerance, the dependent variables are updated and D x is increased to its optimum value.

The embeddedRK5() method source code follows .

 public static int embeddedRK5(ODE ode, double range,                          double dx, double tolerance) {   double maxError;   int i,j,k,m;   //  Define some convenience variables to make   //  the code more readable.   int numEqns = ode.getNumEqns();   double x[] = ode.getX();     double y[][] = ode.getY();   //  Create some local arrays   double dy[][] = new double[6][numEqns];   double dyTotal[] = new double[numEqns];   double ytmp[] = new double[numEqns];   double error[][] =          new double[ODE.MAX_STEPS][numEqns];   //  load the Cash-Karp parameters   double a[] = {0.0, 0.2, 0.3, 0.6, 1.0, 0.875};   double c[] = {37.0/378.0, 0.0, 250.0/621.0,                 125.0/594.0, 0.0, 512.0/1771.0};   double d[] = new double[6];   d[0] = c[0] - 2825.0/27648.0;   d[1] = 0.0;   d[2] = c[2] - 18575.0/48384.0;   d[3] = c[3] - 13525.0/55296.0;   d[4] = c[4] - 277.0/14336.0;   d[5] = c[5] - 0.25;   double b[][] = { {0.0, 0.0, 0.0, 0.0, 0.0},                    {0.2, 0.0, 0.0, 0.0, 0.0},                    {0.075, 0.225, 0.0, 0.0, 0.0},                    {0.3, -0.9, 1.2, 0.0, 0.0},          {-11.0/54.0, 2.5, -70.0/27.0, 35.0/27.0, 0.0},          {1631.0/55296.0, 175.0/512.0, 575.0/13824.0,           44275.0/110592.0, 253.0/4096.0} };   //  Integrate the ODE over the desired range.   //  Stop if you are going to overflow the matrices   i=1;   while( x[i-1] < range && i < ODE.MAX_STEPS-1) {     //  Set up an iteration loop to optimize dx     while (true) {       //  First Runge-Kutta step       ode.getFunction(x[i-1],dy[0],y[i-1]);       for(j=0; j<numEqns; ++j) {         dy[0][j] *= dx;         dyTotal[j] = c[0]*dy[0][j];         error[i][j] = d[0]*dy[0][j];       }       //  Runge-Kutta steps 2-6       for(k=1; k<6; ++k) {         for(j=0; j<numEqns; ++j) {           ytmp[j] = y[i-1][j];           for(m=0; m<k; ++m) {             ytmp[j] += b[k][m]*dy[m][j];           }         }         ode.getFunction(x[i-1]+a[k]*dx,dy[k],ytmp);         for(j=0; j<numEqns; ++j) {           dy[k][j] *= dx;           dyTotal[j] += c[k]*dy[k][j];           error[i][j] += d[k]*dy[k][j];         }       }       //  Compute maximum error       maxError = 0.0;       for(j=0; j<numEqns; ++j) {         maxError =             Math.max(maxError, Math.abs(error[i][j]));       }       //  If the maximum error is greater than the       //  tolerance, decrease delta-x and try again.       //  Otherwise, update the variables and move on to       //  the next point.       if ( maxError > tolerance ) {         dx *= Math.pow(tolerance/maxError,0.2);       } else {         break;       }     }     //  Update the dependent variables     for(j=0; j<numEqns; ++j) {       y[i][j] = y[i-1][j] + dyTotal[j];     }     //  Increment independent variable, reset dx, and     //  move on to the next point. Make sure you don't     //  go past the specified range.     x[i] = x[i-1] + dx;     dx *= Math.pow(tolerance/maxError,0.2);     if ( x[i]+dx > range ) {       dx = range - x[i];     }     //  Go to the next dependent variable location     ++i;   }  //  end of outer while loop   //  Return the number of steps computed   return i; } 

Let's use the embeddedRK5() method to solve the same spring problem that was solved by the rungeKutta4() method in the earlier example. The EmbedSpring.java program is a driver program that creates a SpringODE object with the same initial values as the one from the rungeKutta4() method example. The ODE is solved by calling the embeddedRK5() method. The truncation error tolerance is set to be 1.0e-6.

 import TechJava.MathLib.*; public class EmbedSpring {   public static void main(String args[]) {     //  Create a SpringODE object     double mass = 1.0;     double mu = 1.5;     double k = 20.0;     SpringODE ode = new SpringODE(k, mu, mass);     //  load initial conditions.  The spring is     //  initially stretched 0.1 meters from its     //  equilibrium  position.     double V[] = {-0.2};     ode.setInitialConditions(V);     //  Solve the ODE over the desired range with the     //  specified tolerance and initial step size     double dx = 0.1;     double range = 5.0;     double tolerance = 1.0e-6;     int numSteps =        ODESolver.embeddedRK5(ode, range, dx, tolerance);     //  Print out the results     System.out.println("i   t   dxdt   x");     for(int i=0; i<numSteps; ++i) {       System.out.println(            ""+i+" "+ode.getOneX(i)+            " "+ode.getOneY(i,0)+ " "+ode.getOneY(i,1));     }   } } 

The output of the EmbedSpring program is shown in Figure 20.3. At first glance it seems quite similar to the results generated by the rungeKutta4() method. The output from the embeddedRK5() method tracks the exact solution very closely. If you look at the distribution of points you will notice that the embedded Runge-Kutta algorithm placed more points in regions of high gradients (around the maximum and minimum amplitudes) and fewer points in the smoother regions of the curve.

Figure 20.3. Spring position as a function of time

graphics/20fig03.gif



Technical Java. Applications for Science and Engineering
Technical Java: Applications for Science and Engineering
ISBN: 0131018159
EAN: 2147483647
Year: 2003
Pages: 281
Authors: Grant Palmer

flylib.com © 2008-2017.
If you may any questions please contact us: flylib@qtcs.net