001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.ode.nonstiff;
019    
020    
021    /**
022     * This class implements the 5(4) Higham and Hall integrator for
023     * Ordinary Differential Equations.
024     *
025     * <p>This integrator is an embedded Runge-Kutta integrator
026     * of order 5(4) used in local extrapolation mode (i.e. the solution
027     * is computed using the high order formula) with stepsize control
028     * (and automatic step initialization) and continuous output. This
029     * method uses 7 functions evaluations per step.</p>
030     *
031     * @version $Revision: 810196 $ $Date: 2009-09-01 15:47:46 -0400 (Tue, 01 Sep 2009) $
032     * @since 1.2
033     */
034    
035    public class HighamHall54Integrator extends EmbeddedRungeKuttaIntegrator {
036    
037      /** Integrator method name. */
038      private static final String METHOD_NAME = "Higham-Hall 5(4)";
039    
040      /** Time steps Butcher array. */
041      private static final double[] STATIC_C = {
042        2.0/9.0, 1.0/3.0, 1.0/2.0, 3.0/5.0, 1.0, 1.0
043      };
044    
045      /** Internal weights Butcher array. */
046      private static final double[][] STATIC_A = {
047        {2.0/9.0},
048        {1.0/12.0, 1.0/4.0},
049        {1.0/8.0, 0.0, 3.0/8.0},
050        {91.0/500.0, -27.0/100.0, 78.0/125.0, 8.0/125.0},
051        {-11.0/20.0, 27.0/20.0, 12.0/5.0, -36.0/5.0, 5.0},
052        {1.0/12.0, 0.0, 27.0/32.0, -4.0/3.0, 125.0/96.0, 5.0/48.0}
053      };
054    
055      /** Propagation weights Butcher array. */
056      private static final double[] STATIC_B = {
057        1.0/12.0, 0.0, 27.0/32.0, -4.0/3.0, 125.0/96.0, 5.0/48.0, 0.0
058      };
059    
060      /** Error weights Butcher array. */
061      private static final double[] STATIC_E = {
062        -1.0/20.0, 0.0, 81.0/160.0, -6.0/5.0, 25.0/32.0, 1.0/16.0, -1.0/10.0
063      };
064    
065      /** Simple constructor.
066       * Build a fifth order Higham and Hall integrator with the given step bounds
067       * @param minStep minimal step (must be positive even for backward
068       * integration), the last step can be smaller than this
069       * @param maxStep maximal step (must be positive even for backward
070       * integration)
071       * @param scalAbsoluteTolerance allowed absolute error
072       * @param scalRelativeTolerance allowed relative error
073       */
074      public HighamHall54Integrator(final double minStep, final double maxStep,
075                                    final double scalAbsoluteTolerance,
076                                    final double scalRelativeTolerance) {
077        super(METHOD_NAME, false, STATIC_C, STATIC_A, STATIC_B, new HighamHall54StepInterpolator(),
078              minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
079      }
080    
081      /** Simple constructor.
082       * Build a fifth order Higham and Hall integrator with the given step bounds
083       * @param minStep minimal step (must be positive even for backward
084       * integration), the last step can be smaller than this
085       * @param maxStep maximal step (must be positive even for backward
086       * integration)
087       * @param vecAbsoluteTolerance allowed absolute error
088       * @param vecRelativeTolerance allowed relative error
089       */
090      public HighamHall54Integrator(final double minStep, final double maxStep,
091                                    final double[] vecAbsoluteTolerance,
092                                    final double[] vecRelativeTolerance) {
093        super(METHOD_NAME, false, STATIC_C, STATIC_A, STATIC_B, new HighamHall54StepInterpolator(),
094              minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
095      }
096    
097      /** {@inheritDoc} */
098      @Override
099      public int getOrder() {
100        return 5;
101      }
102    
103      /** {@inheritDoc} */
104      @Override
105      protected double estimateError(final double[][] yDotK,
106                                     final double[] y0, final double[] y1,
107                                     final double h) {
108    
109        double error = 0;
110    
111        for (int j = 0; j < y0.length; ++j) {
112          double errSum = STATIC_E[0] * yDotK[0][j];
113          for (int l = 1; l < STATIC_E.length; ++l) {
114            errSum += STATIC_E[l] * yDotK[l][j];
115          }
116    
117          final double yScale = Math.max(Math.abs(y0[j]), Math.abs(y1[j]));
118          final double tol = (vecAbsoluteTolerance == null) ?
119                             (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
120                             (vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
121          final double ratio  = h * errSum / tol;
122          error += ratio * ratio;
123    
124        }
125    
126        return Math.sqrt(error / y0.length);
127    
128      }
129    
130    }