1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/*
2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more
3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements.  See the NOTICE file distributed with
4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership.
5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0
6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with
7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License.  You may obtain a copy of the License at
8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond *
9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond *      http://www.apache.org/licenses/LICENSE-2.0
10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond *
11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software
12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS,
13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and
15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License.
16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */
17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.ode.nonstiff;
19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.io.IOException;
21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.io.ObjectInput;
22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.io.ObjectOutput;
23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.AbstractIntegrator;
25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** This class represents an interpolator over the last step during an
28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ODE integration for Runge-Kutta and embedded Runge-Kutta integrators.
29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond *
30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see RungeKuttaIntegrator
31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see EmbeddedRungeKuttaIntegrator
32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond *
33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 811827 $ $Date: 2009-09-06 17:32:50 +0200 (dim. 06 sept. 2009) $
34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 1.2
35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */
36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondabstract class RungeKuttaStepInterpolator
38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  extends AbstractStepInterpolator {
39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    /** Slopes at the intermediate points */
41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    protected double[][] yDotK;
42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    /** Reference to the integrator. */
44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    protected AbstractIntegrator integrator;
45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  /** Simple constructor.
47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * This constructor builds an instance that is not usable yet, the
48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * {@link #reinitialize} method should be called before using the
49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * instance in order to initialize the internal arrays. This
50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * constructor is used only in order to delay the initialization in
51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * some cases. The {@link RungeKuttaIntegrator} and {@link
52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * EmbeddedRungeKuttaIntegrator} classes use the prototyping design
53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * pattern to create the step interpolators by cloning an
54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * uninitialized model and latter initializing the copy.
55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   */
56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  protected RungeKuttaStepInterpolator() {
57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    super();
58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    yDotK      = null;
59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    integrator = null;
60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  }
61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  /** Copy constructor.
63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  * <p>The copied interpolator should have been finalized before the
65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  * copy, otherwise the copy will not be able to perform correctly any
66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  * interpolation and will throw a {@link NullPointerException}
67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  * later. Since we don't want this constructor to throw the
68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  * exceptions finalization may involve and since we don't want this
69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  * method to modify the state of the copied interpolator,
70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  * finalization is <strong>not</strong> done automatically, it
71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  * remains under user control.</p>
72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  * <p>The copy is a deep copy: its arrays are separated from the
74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  * original arrays of the instance.</p>
75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  * @param interpolator interpolator to copy from.
77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  */
79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  public RungeKuttaStepInterpolator(final RungeKuttaStepInterpolator interpolator) {
80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    super(interpolator);
82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    if (interpolator.currentState != null) {
84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond      final int dimension = currentState.length;
85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond      yDotK = new double[interpolator.yDotK.length][];
87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond      for (int k = 0; k < interpolator.yDotK.length; ++k) {
88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        yDotK[k] = new double[dimension];
89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        System.arraycopy(interpolator.yDotK[k], 0,
90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond                         yDotK[k], 0, dimension);
91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond      }
92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    } else {
94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond      yDotK = null;
95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    }
96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    // we cannot keep any reference to the equations in the copy
98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    // the interpolator should have been finalized before
99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    integrator = null;
100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  }
102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  /** Reinitialize the instance
104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * <p>Some Runge-Kutta integrators need fewer functions evaluations
105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * than their counterpart step interpolators. So the interpolator
106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * should perform the last evaluations they need by themselves. The
107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * {@link RungeKuttaIntegrator RungeKuttaIntegrator} and {@link
108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * EmbeddedRungeKuttaIntegrator EmbeddedRungeKuttaIntegrator}
109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * abstract classes call this method in order to let the step
110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * interpolator perform the evaluations it needs. These evaluations
111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * will be performed during the call to <code>doFinalize</code> if
112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * any, i.e. only if the step handler either calls the {@link
113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * AbstractStepInterpolator#finalizeStep finalizeStep} method or the
114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * {@link AbstractStepInterpolator#getInterpolatedState
115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * getInterpolatedState} method (for an interpolator which needs a
116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * finalization) or if it clones the step interpolator.</p>
117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * @param rkIntegrator integrator being used
118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * @param y reference to the integrator array holding the state at
119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * the end of the step
120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * @param yDotArray reference to the integrator array holding all the
121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * intermediate slopes
122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   * @param forward integration direction indicator
123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond   */
124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  public void reinitialize(final AbstractIntegrator rkIntegrator,
125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond                           final double[] y, final double[][] yDotArray, final boolean forward) {
126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    reinitialize(y, forward);
127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    this.yDotK = yDotArray;
128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    this.integrator = rkIntegrator;
129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  }
130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  /** {@inheritDoc} */
132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  @Override
133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  public void writeExternal(final ObjectOutput out)
134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    throws IOException {
135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    // save the state of the base class
137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    writeBaseExternal(out);
138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    // save the local attributes
140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    final int n = (currentState == null) ? -1 : currentState.length;
141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    final int kMax = (yDotK == null) ? -1 : yDotK.length;
142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    out.writeInt(kMax);
143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    for (int k = 0; k < kMax; ++k) {
144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond      for (int i = 0; i < n; ++i) {
145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        out.writeDouble(yDotK[k][i]);
146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond      }
147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    }
148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    // we do not save any reference to the equations
150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  }
152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  /** {@inheritDoc} */
154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  @Override
155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  public void readExternal(final ObjectInput in)
156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    throws IOException {
157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    // read the base class
159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    final double t = readBaseExternal(in);
160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    // read the local attributes
162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    final int n = (currentState == null) ? -1 : currentState.length;
163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    final int kMax = in.readInt();
164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    yDotK = (kMax < 0) ? null : new double[kMax][];
165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    for (int k = 0; k < kMax; ++k) {
166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond      yDotK[k] = (n < 0) ? null : new double[n];
167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond      for (int i = 0; i < n; ++i) {
168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        yDotK[k][i] = in.readDouble();
169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond      }
170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    }
171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    integrator = null;
173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    if (currentState != null) {
175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        // we can now set the interpolated time and state
176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        setInterpolatedTime(t);
177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    } else {
178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        interpolatedTime = t;
179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    }
180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond  }
182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond}
184