Home | History | Annotate | Download | only in nonstiff 
      1 /*
      2  * Licensed to the Apache Software Foundation (ASF) under one or more
      3  * contributor license agreements.  See the NOTICE file distributed with
      4  * this work for additional information regarding copyright ownership.
      5  * The ASF licenses this file to You under the Apache License, Version 2.0
      6  * (the "License"); you may not use this file except in compliance with
      7  * the License.  You may obtain a copy of the License at
      8  *
      9  *      http://www.apache.org/licenses/LICENSE-2.0
     10  *
     11  * Unless required by applicable law or agreed to in writing, software
     12  * distributed under the License is distributed on an "AS IS" BASIS,
     13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
     14  * See the License for the specific language governing permissions and
     15  * limitations under the License.
     16  */
     17 
     18 package org.apache.commons.math.ode.nonstiff;
     19 
     20 import org.apache.commons.math.ode.DerivativeException;
     21 import org.apache.commons.math.ode.sampling.StepInterpolator;
     22 
     23 /**
     24  * This class implements a step interpolator for the classical fourth
     25  * order Runge-Kutta integrator.
     26  *
     27  * <p>This interpolator allows to compute dense output inside the last
     28  * step computed. The interpolation equation is consistent with the
     29  * integration scheme :
     30 
     31  * <pre>
     32  *   y(t_n + theta h) = y (t_n + h)
     33  *                    + (1 - theta) (h/6) [ (-4 theta^2 + 5 theta - 1) y'_1
     34  *                                          +(4 theta^2 - 2 theta - 2) (y'_2 + y'_3)
     35  *                                          -(4 theta^2 +   theta + 1) y'_4
     36  *                                        ]
     37  * </pre>
     38  *
     39  * where theta belongs to [0 ; 1] and where y'_1 to y'_4 are the four
     40  * evaluations of the derivatives already computed during the
     41  * step.</p>
     42  *
     43  * @see ClassicalRungeKuttaIntegrator
     44  * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 fvr. 2011) $
     45  * @since 1.2
     46  */
     47 
     48 class ClassicalRungeKuttaStepInterpolator
     49     extends RungeKuttaStepInterpolator {
     50 
     51     /** Serializable version identifier */
     52     private static final long serialVersionUID = -6576285612589783992L;
     53 
     54     /** Simple constructor.
     55      * This constructor builds an instance that is not usable yet, the
     56      * {@link RungeKuttaStepInterpolator#reinitialize} method should be
     57      * called before using the instance in order to initialize the
     58      * internal arrays. This constructor is used only in order to delay
     59      * the initialization in some cases. The {@link RungeKuttaIntegrator}
     60      * class uses the prototyping design pattern to create the step
     61      * interpolators by cloning an uninitialized model and latter initializing
     62      * the copy.
     63      */
     64     public ClassicalRungeKuttaStepInterpolator() {
     65     }
     66 
     67     /** Copy constructor.
     68      * @param interpolator interpolator to copy from. The copy is a deep
     69      * copy: its arrays are separated from the original arrays of the
     70      * instance
     71      */
     72     public ClassicalRungeKuttaStepInterpolator(final ClassicalRungeKuttaStepInterpolator interpolator) {
     73         super(interpolator);
     74     }
     75 
     76     /** {@inheritDoc} */
     77     @Override
     78     protected StepInterpolator doCopy() {
     79         return new ClassicalRungeKuttaStepInterpolator(this);
     80     }
     81 
     82     /** {@inheritDoc} */
     83     @Override
     84     protected void computeInterpolatedStateAndDerivatives(final double theta,
     85                                             final double oneMinusThetaH)
     86         throws DerivativeException {
     87 
     88         final double fourTheta      = 4 * theta;
     89         final double oneMinusTheta  = 1 - theta;
     90         final double oneMinus2Theta = 1 - 2 * theta;
     91         final double s             = oneMinusThetaH / 6.0;
     92         final double coeff1        = s * ((-fourTheta + 5) * theta - 1);
     93         final double coeff23       = s * (( fourTheta - 2) * theta - 2);
     94         final double coeff4        = s * ((-fourTheta - 1) * theta - 1);
     95         final double coeffDot1     = oneMinusTheta * oneMinus2Theta;
     96         final double coeffDot23    = 2 * theta * oneMinusTheta;
     97         final double coeffDot4     = -theta * oneMinus2Theta;
     98         for (int i = 0; i < interpolatedState.length; ++i) {
     99             final double yDot1  = yDotK[0][i];
    100             final double yDot23 = yDotK[1][i] + yDotK[2][i];
    101             final double yDot4  = yDotK[3][i];
    102             interpolatedState[i] =
    103                 currentState[i] + coeff1  * yDot1 + coeff23 * yDot23 + coeff4  * yDot4;
    104             interpolatedDerivatives[i] =
    105                 coeffDot1 * yDot1 + coeffDot23 * yDot23 + coeffDot4 * yDot4;
    106         }
    107 
    108     }
    109 
    110 }
    111