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> + θ 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 * </p>
48 *
49 * where θ 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 }