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