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 import org.apache.commons.math3.util.FastMath;
22 
23 /**
24  * This class represents an interpolator over the last step during an
25  * ODE integration for the 6th order Luther integrator.
26  *
27  * <p>This interpolator computes dense output inside the last
28  * step computed. The interpolation equation is consistent with the
29  * integration scheme.</p>
30  *
31  * @see LutherIntegrator
32  * @since 3.3
33  */
34 
35 class LutherStepInterpolator extends RungeKuttaStepInterpolator {
36 
37     /** Serializable version identifier */
38     private static final long serialVersionUID = 20140416L;
39 
40     /** Square root. */
41     private static final double Q = FastMath.sqrt(21);
42 
43     /** Simple constructor.
44      * This constructor builds an instance that is not usable yet, the
45      * {@link
46      * org.apache.commons.math3.ode.sampling.AbstractStepInterpolator#reinitialize}
47      * method should be called before using the instance in order to
48      * initialize the internal arrays. This constructor is used only
49      * in order to delay the initialization in some cases. The {@link
50      * RungeKuttaIntegrator} class uses the prototyping design pattern
51      * to create the step interpolators by cloning an uninitialized model
52      * and later initializing the copy.
53      */
54     // CHECKSTYLE: stop RedundantModifier
55     // the public modifier here is needed for serialization
LutherStepInterpolator()56     public LutherStepInterpolator() {
57     }
58     // CHECKSTYLE: resume RedundantModifier
59 
60     /** Copy constructor.
61      * @param interpolator interpolator to copy from. The copy is a deep
62      * copy: its arrays are separated from the original arrays of the
63      * instance
64      */
LutherStepInterpolator(final LutherStepInterpolator interpolator)65     LutherStepInterpolator(final LutherStepInterpolator interpolator) {
66         super(interpolator);
67     }
68 
69     /** {@inheritDoc} */
70     @Override
doCopy()71     protected StepInterpolator doCopy() {
72         return new LutherStepInterpolator(this);
73     }
74 
75 
76     /** {@inheritDoc} */
77     @Override
computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH)78     protected void computeInterpolatedStateAndDerivatives(final double theta,
79                                                           final double oneMinusThetaH) {
80 
81         // the coefficients below have been computed by solving the
82         // order conditions from a theorem from Butcher (1963), using
83         // the method explained in Folkmar Bornemann paper "Runge-Kutta
84         // Methods, Trees, and Maple", Center of Mathematical Sciences, Munich
85         // University of Technology, February 9, 2001
86         //<http://wwwzenger.informatik.tu-muenchen.de/selcuk/sjam012101.html>
87 
88         // the method is implemented in the rkcheck tool
89         // <https://www.spaceroots.org/software/rkcheck/index.html>.
90         // Running it for order 5 gives the following order conditions
91         // for an interpolator:
92         // order 1 conditions
93         // \sum_{i=1}^{i=s}\left(b_{i} \right) =1
94         // order 2 conditions
95         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}\right) = \frac{\theta}{2}
96         // order 3 conditions
97         // \sum_{i=2}^{i=s}\left(b_{i} \sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)}\right) = \frac{\theta^{2}}{6}
98         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{2}\right) = \frac{\theta^{2}}{3}
99         // order 4 conditions
100         // \sum_{i=3}^{i=s}\left(b_{i} \sum_{j=2}^{j=i-1}{\left(a_{i,j} \sum_{k=1}^{k=j-1}{\left(a_{j,k} c_{k} \right)} \right)}\right) = \frac{\theta^{3}}{24}
101         // \sum_{i=2}^{i=s}\left(b_{i} \sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j}^{2} \right)}\right) = \frac{\theta^{3}}{12}
102         // \sum_{i=2}^{i=s}\left(b_{i} c_{i}\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)}\right) = \frac{\theta^{3}}{8}
103         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{3}\right) = \frac{\theta^{3}}{4}
104         // order 5 conditions
105         // \sum_{i=4}^{i=s}\left(b_{i} \sum_{j=3}^{j=i-1}{\left(a_{i,j} \sum_{k=2}^{k=j-1}{\left(a_{j,k} \sum_{l=1}^{l=k-1}{\left(a_{k,l} c_{l} \right)} \right)} \right)}\right) = \frac{\theta^{4}}{120}
106         // \sum_{i=3}^{i=s}\left(b_{i} \sum_{j=2}^{j=i-1}{\left(a_{i,j} \sum_{k=1}^{k=j-1}{\left(a_{j,k} c_{k}^{2} \right)} \right)}\right) = \frac{\theta^{4}}{60}
107         // \sum_{i=3}^{i=s}\left(b_{i} \sum_{j=2}^{j=i-1}{\left(a_{i,j} c_{j}\sum_{k=1}^{k=j-1}{\left(a_{j,k} c_{k} \right)} \right)}\right) = \frac{\theta^{4}}{40}
108         // \sum_{i=2}^{i=s}\left(b_{i} \sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j}^{3} \right)}\right) = \frac{\theta^{4}}{20}
109         // \sum_{i=3}^{i=s}\left(b_{i} c_{i}\sum_{j=2}^{j=i-1}{\left(a_{i,j} \sum_{k=1}^{k=j-1}{\left(a_{j,k} c_{k} \right)} \right)}\right) = \frac{\theta^{4}}{30}
110         // \sum_{i=2}^{i=s}\left(b_{i} c_{i}\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j}^{2} \right)}\right) = \frac{\theta^{4}}{15}
111         // \sum_{i=2}^{i=s}\left(b_{i} \left(\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)} \right)^{2}\right) = \frac{\theta^{4}}{20}
112         // \sum_{i=2}^{i=s}\left(b_{i} c_{i}^{2}\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)}\right) = \frac{\theta^{4}}{10}
113         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{4}\right) = \frac{\theta^{4}}{5}
114 
115         // The a_{j,k} and c_{k} are given by the integrator Butcher arrays. What remains to solve
116         // are the b_i for the interpolator. They are found by solving the above equations.
117         // For a given interpolator, some equations are redundant, so in our case when we select
118         // all equations from order 1 to 4, we still don't have enough independent equations
119         // to solve from b_1 to b_7. We need to also select one equation from order 5. Here,
120         // we selected the last equation. It appears this choice implied at least the last 3 equations
121         // are fulfilled, but some of the former ones are not, so the resulting interpolator is order 5.
122         // At the end, we get the b_i as polynomials in theta.
123 
124         final double coeffDot1 =  1 + theta * ( -54            /   5.0 + theta * (   36                   + theta * ( -47                   + theta *   21)));
125         final double coeffDot2 =  0;
126         final double coeffDot3 =      theta * (-208            /  15.0 + theta * (  320            / 3.0  + theta * (-608            /  3.0 + theta *  112)));
127         final double coeffDot4 =      theta * ( 324            /  25.0 + theta * ( -486            / 5.0  + theta * ( 972            /  5.0 + theta * -567           /  5.0)));
128         final double coeffDot5 =      theta * ((833 + 343 * Q) / 150.0 + theta * ((-637 - 357 * Q) / 30.0 + theta * ((392 + 287 * Q) / 15.0 + theta * (-49 - 49 * Q) /  5.0)));
129         final double coeffDot6 =      theta * ((833 - 343 * Q) / 150.0 + theta * ((-637 + 357 * Q) / 30.0 + theta * ((392 - 287 * Q) / 15.0 + theta * (-49 + 49 * Q) /  5.0)));
130         final double coeffDot7 =      theta * (   3            /   5.0 + theta * (   -3                   + theta *     3));
131 
132         if ((previousState != null) && (theta <= 0.5)) {
133 
134             final double coeff1    =  1 + theta * ( -27            /   5.0 + theta * (   12                   + theta * ( -47            /  4.0 + theta *   21           /  5.0)));
135             final double coeff2    =  0;
136             final double coeff3    =      theta * (-104            /  15.0 + theta * (  320            / 9.0  + theta * (-152            /  3.0 + theta *  112           /  5.0)));
137             final double coeff4    =      theta * ( 162            /  25.0 + theta * ( -162            / 5.0  + theta * ( 243            /  5.0 + theta * -567           / 25.0)));
138             final double coeff5    =      theta * ((833 + 343 * Q) / 300.0 + theta * ((-637 - 357 * Q) / 90.0 + theta * ((392 + 287 * Q) / 60.0 + theta * (-49 - 49 * Q) / 25.0)));
139             final double coeff6    =      theta * ((833 - 343 * Q) / 300.0 + theta * ((-637 + 357 * Q) / 90.0 + theta * ((392 - 287 * Q) / 60.0 + theta * (-49 + 49 * Q) / 25.0)));
140             final double coeff7    =      theta * (   3            /  10.0 + theta * (   -1                   + theta * (   3            /  4.0)));
141             for (int i = 0; i < interpolatedState.length; ++i) {
142                 final double yDot1 = yDotK[0][i];
143                 final double yDot2 = yDotK[1][i];
144                 final double yDot3 = yDotK[2][i];
145                 final double yDot4 = yDotK[3][i];
146                 final double yDot5 = yDotK[4][i];
147                 final double yDot6 = yDotK[5][i];
148                 final double yDot7 = yDotK[6][i];
149                 interpolatedState[i] = previousState[i] +
150                         theta * h * (coeff1 * yDot1 + coeff2 * yDot2 + coeff3 * yDot3 +
151                                      coeff4 * yDot4 + coeff5 * yDot5 + coeff6 * yDot6 + coeff7 * yDot7);
152                 interpolatedDerivatives[i] = coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3 * yDot3 +
153                         coeffDot4 * yDot4 + coeffDot5 * yDot5 + coeffDot6 * yDot6 + coeffDot7 * yDot7;
154             }
155         } else {
156 
157             final double coeff1    =  -1 /  20.0 + theta * (  19            /  20.0 + theta * (  -89             /  20.0  + theta * (   151            /  20.0 + theta *  -21           /   5.0)));
158             final double coeff2    =  0;
159             final double coeff3    = -16 /  45.0 + theta * ( -16            /  45.0 + theta * ( -328             /  45.0  + theta * (   424            /  15.0 + theta * -112           /   5.0)));
160             final double coeff4    =               theta * (                          theta * (  162             /  25.0  + theta * (  -648            /  25.0 + theta *  567           /  25.0)));
161             final double coeff5    = -49 / 180.0 + theta * ( -49            / 180.0 + theta * ((2254 + 1029 * Q) / 900.0  + theta * ((-1372 - 847 * Q) / 300.0 + theta * ( 49 + 49 * Q) /  25.0)));
162             final double coeff6    = -49 / 180.0 + theta * ( -49            / 180.0 + theta * ((2254 - 1029 * Q) / 900.0  + theta * ((-1372 + 847 * Q) / 300.0 + theta * ( 49 - 49 * Q) /  25.0)));
163             final double coeff7    =  -1 /  20.0 + theta * (  -1            /  20.0 + theta * (    1             /   4.0  + theta * (    -3            /   4.0)));
164             for (int i = 0; i < interpolatedState.length; ++i) {
165                 final double yDot1 = yDotK[0][i];
166                 final double yDot2 = yDotK[1][i];
167                 final double yDot3 = yDotK[2][i];
168                 final double yDot4 = yDotK[3][i];
169                 final double yDot5 = yDotK[4][i];
170                 final double yDot6 = yDotK[5][i];
171                 final double yDot7 = yDotK[6][i];
172                 interpolatedState[i] = currentState[i] +
173                         oneMinusThetaH * (coeff1 * yDot1 + coeff2 * yDot2 + coeff3 * yDot3 +
174                                           coeff4 * yDot4 + coeff5 * yDot5 + coeff6 * yDot6 + coeff7 * yDot7);
175                 interpolatedDerivatives[i] = coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3 * yDot3 +
176                         coeffDot4 * yDot4 + coeffDot5 * yDot5 + coeffDot6 * yDot6 + coeffDot7 * yDot7;
177             }
178         }
179 
180     }
181 
182 }
183