View Javadoc
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.math3.ode.nonstiff;
19  
20  import org.apache.commons.math3.ode.sampling.StepInterpolator;
21  
22  /**
23   * This class implements a step interpolator for the 3/8 fourth
24   * order Runge-Kutta integrator.
25   *
26   * <p>This interpolator allows to compute dense output inside the last
27   * step computed. The interpolation equation is consistent with the
28   * integration scheme :
29   * <ul>
30   *   <li>Using reference point at step start:<br>
31   *     y(t<sub>n</sub> + &theta; h) = y (t<sub>n</sub>)
32   *                      + &theta; (h/8) [ (8 - 15 &theta; +  8 &theta;<sup>2</sup>) y'<sub>1</sub>
33   *                                     +  3 * (15 &theta; - 12 &theta;<sup>2</sup>) y'<sub>2</sub>
34   *                                     +        3 &theta;                           y'<sub>3</sub>
35   *                                     +      (-3 &theta; +  4 &theta;<sup>2</sup>) y'<sub>4</sub>
36   *                                    ]
37   *   </li>
38   *   <li>Using reference point at step end:<br>
39   *     y(t<sub>n</sub> + &theta; h) = y (t<sub>n</sub> + h)
40   *                      - (1 - &theta;) (h/8) [(1 - 7 &theta; + 8 &theta;<sup>2</sup>) y'<sub>1</sub>
41   *                                         + 3 (1 +   &theta; - 4 &theta;<sup>2</sup>) y'<sub>2</sub>
42   *                                         + 3 (1 +   &theta;)                         y'<sub>3</sub>
43   *                                         +   (1 +   &theta; + 4 &theta;<sup>2</sup>) y'<sub>4</sub>
44   *                                          ]
45   *   </li>
46   * </ul>
47   * </p>
48   *
49   * where &theta; belongs to [0 ; 1] and where y'<sub>1</sub> to y'<sub>4</sub> are the four
50   * evaluations of the derivatives already computed during the
51   * step.</p>
52   *
53   * @see ThreeEighthesIntegrator
54   * @version $Id: ThreeEighthesStepInterpolator.java 1416643 2012-12-03 19:37:14Z tn $
55   * @since 1.2
56   */
57  
58  class ThreeEighthesStepInterpolator
59    extends RungeKuttaStepInterpolator {
60  
61    /** Serializable version identifier */
62    private static final long serialVersionUID = 20111120L;
63  
64    /** Simple constructor.
65     * This constructor builds an instance that is not usable yet, the
66     * {@link
67     * org.apache.commons.math3.ode.sampling.AbstractStepInterpolator#reinitialize}
68     * method should be called before using the instance in order to
69     * initialize the internal arrays. This constructor is used only
70     * in order to delay the initialization in some cases. The {@link
71     * RungeKuttaIntegrator} class uses the prototyping design pattern
72     * to create the step interpolators by cloning an uninitialized model
73     * and later initializing the copy.
74     */
75    public ThreeEighthesStepInterpolator() {
76    }
77  
78    /** Copy constructor.
79     * @param interpolator interpolator to copy from. The copy is a deep
80     * copy: its arrays are separated from the original arrays of the
81     * instance
82     */
83    public ThreeEighthesStepInterpolator(final ThreeEighthesStepInterpolator interpolator) {
84      super(interpolator);
85    }
86  
87    /** {@inheritDoc} */
88    @Override
89    protected StepInterpolator doCopy() {
90      return new ThreeEighthesStepInterpolator(this);
91    }
92  
93  
94    /** {@inheritDoc} */
95    @Override
96    protected void computeInterpolatedStateAndDerivatives(final double theta,
97                                            final double oneMinusThetaH) {
98  
99        final double coeffDot3  = 0.75 * theta;
100       final double coeffDot1  = coeffDot3 * (4 * theta - 5) + 1;
101       final double coeffDot2  = coeffDot3 * (5 - 6 * theta);
102       final double coeffDot4  = coeffDot3 * (2 * theta - 1);
103 
104       if ((previousState != null) && (theta <= 0.5)) {
105           final double s          = theta * h / 8.0;
106           final double fourTheta2 = 4 * theta * theta;
107           final double coeff1     = s * (8 - 15 * theta + 2 * fourTheta2);
108           final double coeff2     = 3 * s * (5 * theta - fourTheta2);
109           final double coeff3     = 3 * s * theta;
110           final double coeff4     = s * (-3 * theta + fourTheta2);
111           for (int i = 0; i < interpolatedState.length; ++i) {
112               final double yDot1 = yDotK[0][i];
113               final double yDot2 = yDotK[1][i];
114               final double yDot3 = yDotK[2][i];
115               final double yDot4 = yDotK[3][i];
116               interpolatedState[i] =
117                       previousState[i] + coeff1 * yDot1 + coeff2 * yDot2 + coeff3 * yDot3 + coeff4 * yDot4;
118               interpolatedDerivatives[i] =
119                       coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3 * yDot3 + coeffDot4 * yDot4;
120 
121           }
122       } else {
123           final double s          = oneMinusThetaH / 8.0;
124           final double fourTheta2 = 4 * theta * theta;
125           final double coeff1     = s * (1 - 7 * theta + 2 * fourTheta2);
126           final double coeff2     = 3 * s * (1 + theta - fourTheta2);
127           final double coeff3     = 3 * s * (1 + theta);
128           final double coeff4     = s * (1 + theta + fourTheta2);
129           for (int i = 0; i < interpolatedState.length; ++i) {
130               final double yDot1 = yDotK[0][i];
131               final double yDot2 = yDotK[1][i];
132               final double yDot3 = yDotK[2][i];
133               final double yDot4 = yDotK[3][i];
134               interpolatedState[i] =
135                       currentState[i] - coeff1 * yDot1 - coeff2 * yDot2 - coeff3 * yDot3 - coeff4 * yDot4;
136               interpolatedDerivatives[i] =
137                       coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3 * yDot3 + coeffDot4 * yDot4;
138 
139           }
140       }
141 
142   }
143 
144 }