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   * @since 1.2
55   */
56  
57  class ThreeEighthesStepInterpolator
58    extends RungeKuttaStepInterpolator {
59  
60    /** Serializable version identifier */
61    private static final long serialVersionUID = 20111120L;
62  
63    /** Simple constructor.
64     * This constructor builds an instance that is not usable yet, the
65     * {@link
66     * org.apache.commons.math3.ode.sampling.AbstractStepInterpolator#reinitialize}
67     * method should be called before using the instance in order to
68     * initialize the internal arrays. This constructor is used only
69     * in order to delay the initialization in some cases. The {@link
70     * RungeKuttaIntegrator} class uses the prototyping design pattern
71     * to create the step interpolators by cloning an uninitialized model
72     * and later initializing the copy.
73     */
74    public ThreeEighthesStepInterpolator() {
75    }
76  
77    /** Copy constructor.
78     * @param interpolator interpolator to copy from. The copy is a deep
79     * copy: its arrays are separated from the original arrays of the
80     * instance
81     */
82    public ThreeEighthesStepInterpolator(final ThreeEighthesStepInterpolator interpolator) {
83      super(interpolator);
84    }
85  
86    /** {@inheritDoc} */
87    @Override
88    protected StepInterpolator doCopy() {
89      return new ThreeEighthesStepInterpolator(this);
90    }
91  
92  
93    /** {@inheritDoc} */
94    @Override
95    protected void computeInterpolatedStateAndDerivatives(final double theta,
96                                            final double oneMinusThetaH) {
97  
98        final double coeffDot3  = 0.75 * theta;
99        final double coeffDot1  = coeffDot3 * (4 * theta - 5) + 1;
100       final double coeffDot2  = coeffDot3 * (5 - 6 * theta);
101       final double coeffDot4  = coeffDot3 * (2 * theta - 1);
102 
103       if ((previousState != null) && (theta <= 0.5)) {
104           final double s          = theta * h / 8.0;
105           final double fourTheta2 = 4 * theta * theta;
106           final double coeff1     = s * (8 - 15 * theta + 2 * fourTheta2);
107           final double coeff2     = 3 * s * (5 * theta - fourTheta2);
108           final double coeff3     = 3 * s * theta;
109           final double coeff4     = s * (-3 * theta + fourTheta2);
110           for (int i = 0; i < interpolatedState.length; ++i) {
111               final double yDot1 = yDotK[0][i];
112               final double yDot2 = yDotK[1][i];
113               final double yDot3 = yDotK[2][i];
114               final double yDot4 = yDotK[3][i];
115               interpolatedState[i] =
116                       previousState[i] + coeff1 * yDot1 + coeff2 * yDot2 + coeff3 * yDot3 + coeff4 * yDot4;
117               interpolatedDerivatives[i] =
118                       coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3 * yDot3 + coeffDot4 * yDot4;
119 
120           }
121       } else {
122           final double s          = oneMinusThetaH / 8.0;
123           final double fourTheta2 = 4 * theta * theta;
124           final double coeff1     = s * (1 - 7 * theta + 2 * fourTheta2);
125           final double coeff2     = 3 * s * (1 + theta - fourTheta2);
126           final double coeff3     = 3 * s * (1 + theta);
127           final double coeff4     = s * (1 + theta + fourTheta2);
128           for (int i = 0; i < interpolatedState.length; ++i) {
129               final double yDot1 = yDotK[0][i];
130               final double yDot2 = yDotK[1][i];
131               final double yDot3 = yDotK[2][i];
132               final double yDot4 = yDotK[3][i];
133               interpolatedState[i] =
134                       currentState[i] - coeff1 * yDot1 - coeff2 * yDot2 - coeff3 * yDot3 - coeff4 * yDot4;
135               interpolatedDerivatives[i] =
136                       coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3 * yDot3 + coeffDot4 * yDot4;
137 
138           }
139       }
140 
141   }
142 
143 }