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.core.Field;
21 import org.apache.commons.math4.legacy.core.RealFieldElement;
22 import org.apache.commons.math4.legacy.ode.FieldEquationsMapper;
23 import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
24 import org.apache.commons.math4.legacy.core.MathArrays;
25
26
27 /**
28 * This class implements the Luther sixth order Runge-Kutta
29 * integrator for Ordinary Differential Equations.
30
31 * <p>
32 * This method is described in H. A. Luther 1968 paper <a
33 * href="http://www.ams.org/journals/mcom/1968-22-102/S0025-5718-68-99876-1/S0025-5718-68-99876-1.pdf">
34 * An explicit Sixth-Order Runge-Kutta Formula</a>.
35 * </p>
36
37 * <p>This method is an explicit Runge-Kutta method, its Butcher-array
38 * is the following one :
39 * <pre>
40 * 0 | 0 0 0 0 0 0
41 * 1 | 1 0 0 0 0 0
42 * 1/2 | 3/8 1/8 0 0 0 0
43 * 2/3 | 8/27 2/27 8/27 0 0 0
44 * (7-q)/14 | ( -21 + 9q)/392 ( -56 + 8q)/392 ( 336 - 48q)/392 ( -63 + 3q)/392 0 0
45 * (7+q)/14 | (-1155 - 255q)/1960 ( -280 - 40q)/1960 ( 0 - 320q)/1960 ( 63 + 363q)/1960 ( 2352 + 392q)/1960 0
46 * 1 | ( 330 + 105q)/180 ( 120 + 0q)/180 ( -200 + 280q)/180 ( 126 - 189q)/180 ( -686 - 126q)/180 ( 490 - 70q)/180
47 * |--------------------------------------------------------------------------------------------------------------------------------------------------
48 * | 1/20 0 16/45 0 49/180 49/180 1/20
49 * </pre>
50 * where q = √21
51 *
52 * @see EulerFieldIntegrator
53 * @see ClassicalRungeKuttaFieldIntegrator
54 * @see GillFieldIntegrator
55 * @see MidpointFieldIntegrator
56 * @see ThreeEighthesFieldIntegrator
57 * @param <T> the type of the field elements
58 * @since 3.6
59 */
60
61 public class LutherFieldIntegrator<T extends RealFieldElement<T>>
62 extends RungeKuttaFieldIntegrator<T> {
63
64 /** Simple constructor.
65 * Build a fourth-order Luther integrator with the given step.
66 * @param field field to which the time and state vector elements belong
67 * @param step integration step
68 */
69 public LutherFieldIntegrator(final Field<T> field, final T step) {
70 super(field, "Luther", step);
71 }
72
73 /** {@inheritDoc} */
74 @Override
75 public T[] getC() {
76 final T q = getField().getZero().add(21).sqrt();
77 final T[] c = MathArrays.buildArray(getField(), 6);
78 c[0] = getField().getOne();
79 c[1] = fraction(1, 2);
80 c[2] = fraction(2, 3);
81 c[3] = q.subtract(7).divide(-14);
82 c[4] = q.add(7).divide(14);
83 c[5] = getField().getOne();
84 return c;
85 }
86
87 /** {@inheritDoc} */
88 @Override
89 public T[][] getA() {
90 final T q = getField().getZero().add(21).sqrt();
91 final T[][] a = MathArrays.buildArray(getField(), 6, -1);
92 for (int i = 0; i < a.length; ++i) {
93 a[i] = MathArrays.buildArray(getField(), i + 1);
94 }
95 a[0][0] = getField().getOne();
96 a[1][0] = fraction(3, 8);
97 a[1][1] = fraction(1, 8);
98 a[2][0] = fraction(8, 27);
99 a[2][1] = fraction(2, 27);
100 a[2][2] = a[2][0];
101 a[3][0] = q.multiply( 9).add( -21).divide( 392);
102 a[3][1] = q.multiply( 8).add( -56).divide( 392);
103 a[3][2] = q.multiply( -48).add( 336).divide( 392);
104 a[3][3] = q.multiply( 3).add( -63).divide( 392);
105 a[4][0] = q.multiply(-255).add(-1155).divide(1960);
106 a[4][1] = q.multiply( -40).add( -280).divide(1960);
107 a[4][2] = q.multiply(-320) .divide(1960);
108 a[4][3] = q.multiply( 363).add( 63).divide(1960);
109 a[4][4] = q.multiply( 392).add( 2352).divide(1960);
110 a[5][0] = q.multiply( 105).add( 330).divide( 180);
111 a[5][1] = fraction(2, 3);
112 a[5][2] = q.multiply( 280).add( -200).divide( 180);
113 a[5][3] = q.multiply(-189).add( 126).divide( 180);
114 a[5][4] = q.multiply(-126).add( -686).divide( 180);
115 a[5][5] = q.multiply( -70).add( 490).divide( 180);
116 return a;
117 }
118
119 /** {@inheritDoc} */
120 @Override
121 public T[] getB() {
122
123 final T[] b = MathArrays.buildArray(getField(), 7);
124 b[0] = fraction( 1, 20);
125 b[1] = getField().getZero();
126 b[2] = fraction(16, 45);
127 b[3] = getField().getZero();
128 b[4] = fraction(49, 180);
129 b[5] = b[4];
130 b[6] = b[0];
131
132 return b;
133 }
134
135 /** {@inheritDoc} */
136 @Override
137 protected LutherFieldStepInterpolator<T>
138 createInterpolator(final boolean forward, T[][] yDotK,
139 final FieldODEStateAndDerivative<T> globalPreviousState,
140 final FieldODEStateAndDerivative<T> globalCurrentState,
141 final FieldEquationsMapper<T> mapper) {
142 return new LutherFieldStepInterpolator<>(getField(), forward, yDotK,
143 globalPreviousState, globalCurrentState,
144 globalPreviousState, globalCurrentState,
145 mapper);
146 }
147 }