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    // CHECKSTYLE: stop RedundantModifier
75    // the public modifier here is needed for serialization
76    public ThreeEighthesStepInterpolator() {
77    }
78    // CHECKSTYLE: resume RedundantModifier
79  
80    /** Copy constructor.
81     * @param interpolator interpolator to copy from. The copy is a deep
82     * copy: its arrays are separated from the original arrays of the
83     * instance
84     */
85    ThreeEighthesStepInterpolator(final ThreeEighthesStepInterpolator interpolator) {
86      super(interpolator);
87    }
88  
89    /** {@inheritDoc} */
90    @Override
91    protected StepInterpolator doCopy() {
92      return new ThreeEighthesStepInterpolator(this);
93    }
94  
95  
96    /** {@inheritDoc} */
97    @Override
98    protected void computeInterpolatedStateAndDerivatives(final double theta,
99                                            final double oneMinusThetaH) {
100 
101       final double coeffDot3  = 0.75 * theta;
102       final double coeffDot1  = coeffDot3 * (4 * theta - 5) + 1;
103       final double coeffDot2  = coeffDot3 * (5 - 6 * theta);
104       final double coeffDot4  = coeffDot3 * (2 * theta - 1);
105 
106       if ((previousState != null) && (theta <= 0.5)) {
107           final double s          = theta * h / 8.0;
108           final double fourTheta2 = 4 * theta * theta;
109           final double coeff1     = s * (8 - 15 * theta + 2 * fourTheta2);
110           final double coeff2     = 3 * s * (5 * theta - fourTheta2);
111           final double coeff3     = 3 * s * theta;
112           final double coeff4     = s * (-3 * theta + fourTheta2);
113           for (int i = 0; i < interpolatedState.length; ++i) {
114               final double yDot1 = yDotK[0][i];
115               final double yDot2 = yDotK[1][i];
116               final double yDot3 = yDotK[2][i];
117               final double yDot4 = yDotK[3][i];
118               interpolatedState[i] =
119                       previousState[i] + coeff1 * yDot1 + coeff2 * yDot2 + coeff3 * yDot3 + coeff4 * yDot4;
120               interpolatedDerivatives[i] =
121                       coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3 * yDot3 + coeffDot4 * yDot4;
122 
123           }
124       } else {
125           final double s          = oneMinusThetaH / 8.0;
126           final double fourTheta2 = 4 * theta * theta;
127           final double coeff1     = s * (1 - 7 * theta + 2 * fourTheta2);
128           final double coeff2     = 3 * s * (1 + theta - fourTheta2);
129           final double coeff3     = 3 * s * (1 + theta);
130           final double coeff4     = s * (1 + theta + fourTheta2);
131           for (int i = 0; i < interpolatedState.length; ++i) {
132               final double yDot1 = yDotK[0][i];
133               final double yDot2 = yDotK[1][i];
134               final double yDot3 = yDotK[2][i];
135               final double yDot4 = yDotK[3][i];
136               interpolatedState[i] =
137                       currentState[i] - coeff1 * yDot1 - coeff2 * yDot2 - coeff3 * yDot3 - coeff4 * yDot4;
138               interpolatedDerivatives[i] =
139                       coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3 * yDot3 + coeffDot4 * yDot4;
140 
141           }
142       }
143 
144   }
145 
146 }