View Javadoc
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 Gill fourth order Runge-Kutta
29   * integrator for Ordinary Differential Equations .
30  
31   * <p>This method is an explicit Runge-Kutta method, its Butcher-array
32   * is the following one :
33   * <pre>
34   *    0  |    0        0       0      0
35   *   1/2 |   1/2       0       0      0
36   *   1/2 | (q-1)/2  (2-q)/2    0      0
37   *    1  |    0       -q/2  (2+q)/2   0
38   *       |-------------------------------
39   *       |   1/6    (2-q)/6 (2+q)/6  1/6
40   * </pre>
41   * where q = sqrt(2)
42   *
43   * @see EulerFieldIntegrator
44   * @see ClassicalRungeKuttaFieldIntegrator
45   * @see MidpointFieldIntegrator
46   * @see ThreeEighthesFieldIntegrator
47   * @see LutherFieldIntegrator
48   * @param <T> the type of the field elements
49   * @since 3.6
50   */
51  
52  public class GillFieldIntegrator<T extends RealFieldElement<T>>
53      extends RungeKuttaFieldIntegrator<T> {
54  
55      /** Simple constructor.
56       * Build a fourth-order Gill integrator with the given step.
57       * @param field field to which the time and state vector elements belong
58       * @param step integration step
59       */
60      public GillFieldIntegrator(final Field<T> field, final T step) {
61          super(field, "Gill", step);
62      }
63  
64      /** {@inheritDoc} */
65      @Override
66      public T[] getC() {
67          final T[] c = MathArrays.buildArray(getField(), 3);
68          c[0] = fraction(1, 2);
69          c[1] = c[0];
70          c[2] = getField().getOne();
71          return c;
72      }
73  
74      /** {@inheritDoc} */
75      @Override
76      public T[][] getA() {
77  
78          final T two     = getField().getZero().add(2);
79          final T sqrtTwo = two.sqrt();
80  
81          final T[][] a = MathArrays.buildArray(getField(), 3, -1);
82          for (int i = 0; i < a.length; ++i) {
83              a[i] = MathArrays.buildArray(getField(), i + 1);
84          }
85          a[0][0] = fraction(1, 2);
86          a[1][0] = sqrtTwo.subtract(1).multiply(0.5);
87          a[1][1] = sqrtTwo.subtract(2).multiply(-0.5);
88          a[2][0] = getField().getZero();
89          a[2][1] = sqrtTwo.multiply(-0.5);
90          a[2][2] = sqrtTwo.add(2).multiply(0.5);
91          return a;
92      }
93  
94      /** {@inheritDoc} */
95      @Override
96      public T[] getB() {
97  
98          final T two     = getField().getZero().add(2);
99          final T sqrtTwo = two.sqrt();
100 
101         final T[] b = MathArrays.buildArray(getField(), 4);
102         b[0] = fraction(1, 6);
103         b[1] = sqrtTwo.subtract(2).divide(-6);
104         b[2] = sqrtTwo.add(2).divide(6);
105         b[3] = b[0];
106 
107         return b;
108     }
109 
110     /** {@inheritDoc} */
111     @Override
112     protected GillFieldStepInterpolator<T>
113         createInterpolator(final boolean forward, T[][] yDotK,
114                            final FieldODEStateAndDerivative<T> globalPreviousState,
115                            final FieldODEStateAndDerivative<T> globalCurrentState,
116                            final FieldEquationsMapper<T> mapper) {
117         return new GillFieldStepInterpolator<>(getField(), forward, yDotK,
118                                                 globalPreviousState, globalCurrentState,
119                                                 globalPreviousState, globalCurrentState,
120                                                 mapper);
121     }
122 }