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.math4.legacy.ode.nonstiff;
19
20 import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
21
22 /**
23 * This class implements a step interpolator for the classical 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> + θ h) = y (t<sub>n</sub>)
32 * + θ (h/6) [ (6 - 9 θ + 4 θ<sup>2</sup>) y'<sub>1</sub>
33 * + ( 6 θ - 4 θ<sup>2</sup>) (y'<sub>2</sub> + y'<sub>3</sub>)
34 * + ( -3 θ + 4 θ<sup>2</sup>) y'<sub>4</sub>
35 * ]
36 * </li>
37 * <li>Using reference point at step end:<br>
38 * y(t<sub>n</sub> + θ h) = y (t<sub>n</sub> + h)
39 * + (1 - θ) (h/6) [ (-4 θ^2 + 5 θ - 1) y'<sub>1</sub>
40 * +(4 θ^2 - 2 θ - 2) (y'<sub>2</sub> + y'<sub>3</sub>)
41 * -(4 θ^2 + θ + 1) y'<sub>4</sub>
42 * ]
43 * </li>
44 * </ul>
45 *
46 * where θ belongs to [0 ; 1] and where y'<sub>1</sub> to y'<sub>4</sub> are the four
47 * evaluations of the derivatives already computed during the
48 * step.
49 *
50 * @see ClassicalRungeKuttaIntegrator
51 * @since 1.2
52 */
53
54 class ClassicalRungeKuttaStepInterpolator
55 extends RungeKuttaStepInterpolator {
56
57 /** Serializable version identifier. */
58 private static final long serialVersionUID = 20111120L;
59
60 /** Simple constructor.
61 * This constructor builds an instance that is not usable yet, the
62 * {@link RungeKuttaStepInterpolator#reinitialize} method should be
63 * called before using the instance in order to initialize the
64 * internal arrays. This constructor is used only in order to delay
65 * the initialization in some cases. The {@link RungeKuttaIntegrator}
66 * class uses the prototyping design pattern to create the step
67 * interpolators by cloning an uninitialized model and latter initializing
68 * the copy.
69 */
70 // CHECKSTYLE: stop RedundantModifier
71 // the public modifier here is needed for serialization
72 public ClassicalRungeKuttaStepInterpolator() {
73 }
74 // CHECKSTYLE: resume RedundantModifier
75
76 /** Copy constructor.
77 * @param interpolator interpolator to copy from. The copy is a deep
78 * copy: its arrays are separated from the original arrays of the
79 * instance
80 */
81 ClassicalRungeKuttaStepInterpolator(final ClassicalRungeKuttaStepInterpolator interpolator) {
82 super(interpolator);
83 }
84
85 /** {@inheritDoc} */
86 @Override
87 protected StepInterpolator doCopy() {
88 return new ClassicalRungeKuttaStepInterpolator(this);
89 }
90
91 /** {@inheritDoc} */
92 @Override
93 protected void computeInterpolatedStateAndDerivatives(final double theta,
94 final double oneMinusThetaH) {
95
96 final double oneMinusTheta = 1 - theta;
97 final double oneMinus2Theta = 1 - 2 * theta;
98 final double coeffDot1 = oneMinusTheta * oneMinus2Theta;
99 final double coeffDot23 = 2 * theta * oneMinusTheta;
100 final double coeffDot4 = -theta * oneMinus2Theta;
101 if (previousState != null && theta <= 0.5) {
102 final double fourTheta2 = 4 * theta * theta;
103 final double s = theta * h / 6.0;
104 final double coeff1 = s * ( 6 - 9 * theta + fourTheta2);
105 final double coeff23 = s * ( 6 * theta - fourTheta2);
106 final double coeff4 = s * (-3 * theta + fourTheta2);
107 for (int i = 0; i < interpolatedState.length; ++i) {
108 final double yDot1 = yDotK[0][i];
109 final double yDot23 = yDotK[1][i] + yDotK[2][i];
110 final double yDot4 = yDotK[3][i];
111 interpolatedState[i] =
112 previousState[i] + coeff1 * yDot1 + coeff23 * yDot23 + coeff4 * yDot4;
113 interpolatedDerivatives[i] =
114 coeffDot1 * yDot1 + coeffDot23 * yDot23 + coeffDot4 * yDot4;
115 }
116 } else {
117 final double fourTheta = 4 * theta;
118 final double s = oneMinusThetaH / 6.0;
119 final double coeff1 = s * ((-fourTheta + 5) * theta - 1);
120 final double coeff23 = s * (( fourTheta - 2) * theta - 2);
121 final double coeff4 = s * ((-fourTheta - 1) * theta - 1);
122 for (int i = 0; i < interpolatedState.length; ++i) {
123 final double yDot1 = yDotK[0][i];
124 final double yDot23 = yDotK[1][i] + yDotK[2][i];
125 final double yDot4 = yDotK[3][i];
126 interpolatedState[i] =
127 currentState[i] + coeff1 * yDot1 + coeff23 * yDot23 + coeff4 * yDot4;
128 interpolatedDerivatives[i] =
129 coeffDot1 * yDot1 + coeffDot23 * yDot23 + coeffDot4 * yDot4;
130 }
131 }
132 }
133 }