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.math.ode;
19  
20  /**
21   * This class implements a step interpolator for the Gill fourth
22   * order Runge-Kutta integrator.
23   *
24   * <p>This interpolator allows to compute dense output inside the last
25   * step computed. The interpolation equation is consistent with the
26   * integration scheme :
27   *
28   * <pre>
29   *   y(t_n + theta h) = y (t_n + h)
30   *                    - (1 - theta) (h/6) [ (1 - theta) (1 - 4 theta) y'_1
31   *                                        + (1 - theta) (1 + 2 theta) ((2-q) y'_2 + (2+q) y'_3)
32   *                                        + (1 + theta + 4 theta^2) y'_4
33   *                                        ]
34   * </pre>
35   * where theta belongs to [0 ; 1], q = sqrt(2) and where y'_1 to y'_4
36   * are the four evaluations of the derivatives already computed during
37   * the step.</p>
38   *
39   * @see GillIntegrator
40   * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $
41   * @since 1.2
42   */
43  
44  class GillStepInterpolator
45    extends RungeKuttaStepInterpolator {
46      
47    /** Simple constructor.
48     * This constructor builds an instance that is not usable yet, the
49     * {@link AbstractStepInterpolator#reinitialize} method should be called
50     * before using the instance in order to initialize the internal arrays. This
51     * constructor is used only in order to delay the initialization in
52     * some cases. The {@link RungeKuttaIntegrator} class uses the
53     * prototyping design pattern to create the step interpolators by
54     * cloning an uninitialized model and latter initializing the copy.
55     */
56    public GillStepInterpolator() {
57    }
58  
59    /** Copy constructor.
60     * @param interpolator interpolator to copy from. The copy is a deep
61     * copy: its arrays are separated from the original arrays of the
62     * instance
63     */
64    public GillStepInterpolator(GillStepInterpolator interpolator) {
65      super(interpolator);
66    }
67  
68    /** Really copy the finalized instance.
69     * @return a copy of the finalized instance
70     */
71    protected StepInterpolator doCopy() {
72      return new GillStepInterpolator(this);
73    }
74  
75  
76    /** Compute the state at the interpolated time.
77     * This is the main processing method that should be implemented by
78     * the derived classes to perform the interpolation.
79     * @param theta normalized interpolation abscissa within the step
80     * (theta is zero at the previous time step and one at the current time step)
81     * @param oneMinusThetaH time gap between the interpolated time and
82     * the current time
83     * @throws DerivativeException this exception is propagated to the caller if the
84     * underlying user function triggers one
85     */
86    protected void computeInterpolatedState(double theta,
87                                            double oneMinusThetaH)
88      throws DerivativeException {
89  
90      double fourTheta = 4 * theta;
91      double s         = oneMinusThetaH / 6.0;
92      double soMt      = s * (1 - theta);
93      double c23       = soMt * (1 + 2 * theta);
94      double coeff1    = soMt * (1 - fourTheta);
95      double coeff2    = c23  * tMq;
96      double coeff3    = c23  * tPq;
97      double coeff4    = s * (1 + theta * (1 + fourTheta));
98  
99      for (int i = 0; i < interpolatedState.length; ++i) {
100       interpolatedState[i] = currentState[i] -
101                              coeff1 * yDotK[0][i] - coeff2 * yDotK[1][i] -
102                              coeff3 * yDotK[2][i] - coeff4 * yDotK[3][i];
103      }
104 
105   }
106 
107   /** First Gill coefficient. */
108   private static final double tMq = 2 - Math.sqrt(2.0);
109 
110   /** Second Gill coefficient. */
111   private static final double tPq = 2 + Math.sqrt(2.0);
112 
113   /** Serializable version identifier */
114   private static final long serialVersionUID = -107804074496313322L;
115 
116 }