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 18package org.apache.commons.math.ode.nonstiff; 19 20import org.apache.commons.math.ode.DerivativeException; 21import org.apache.commons.math.ode.sampling.StepInterpolator; 22 23/** 24 * This class implements a step interpolator for the classical fourth 25 * order Runge-Kutta integrator. 26 * 27 * <p>This interpolator allows to compute dense output inside the last 28 * step computed. The interpolation equation is consistent with the 29 * integration scheme : 30 31 * <pre> 32 * y(t_n + theta h) = y (t_n + h) 33 * + (1 - theta) (h/6) [ (-4 theta^2 + 5 theta - 1) y'_1 34 * +(4 theta^2 - 2 theta - 2) (y'_2 + y'_3) 35 * -(4 theta^2 + theta + 1) y'_4 36 * ] 37 * </pre> 38 * 39 * where theta belongs to [0 ; 1] and where y'_1 to y'_4 are the four 40 * evaluations of the derivatives already computed during the 41 * step.</p> 42 * 43 * @see ClassicalRungeKuttaIntegrator 44 * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $ 45 * @since 1.2 46 */ 47 48class ClassicalRungeKuttaStepInterpolator 49 extends RungeKuttaStepInterpolator { 50 51 /** Serializable version identifier */ 52 private static final long serialVersionUID = -6576285612589783992L; 53 54 /** Simple constructor. 55 * This constructor builds an instance that is not usable yet, the 56 * {@link RungeKuttaStepInterpolator#reinitialize} method should be 57 * called before using the instance in order to initialize the 58 * internal arrays. This constructor is used only in order to delay 59 * the initialization in some cases. The {@link RungeKuttaIntegrator} 60 * class uses the prototyping design pattern to create the step 61 * interpolators by cloning an uninitialized model and latter initializing 62 * the copy. 63 */ 64 public ClassicalRungeKuttaStepInterpolator() { 65 } 66 67 /** Copy constructor. 68 * @param interpolator interpolator to copy from. The copy is a deep 69 * copy: its arrays are separated from the original arrays of the 70 * instance 71 */ 72 public ClassicalRungeKuttaStepInterpolator(final ClassicalRungeKuttaStepInterpolator interpolator) { 73 super(interpolator); 74 } 75 76 /** {@inheritDoc} */ 77 @Override 78 protected StepInterpolator doCopy() { 79 return new ClassicalRungeKuttaStepInterpolator(this); 80 } 81 82 /** {@inheritDoc} */ 83 @Override 84 protected void computeInterpolatedStateAndDerivatives(final double theta, 85 final double oneMinusThetaH) 86 throws DerivativeException { 87 88 final double fourTheta = 4 * theta; 89 final double oneMinusTheta = 1 - theta; 90 final double oneMinus2Theta = 1 - 2 * theta; 91 final double s = oneMinusThetaH / 6.0; 92 final double coeff1 = s * ((-fourTheta + 5) * theta - 1); 93 final double coeff23 = s * (( fourTheta - 2) * theta - 2); 94 final double coeff4 = s * ((-fourTheta - 1) * theta - 1); 95 final double coeffDot1 = oneMinusTheta * oneMinus2Theta; 96 final double coeffDot23 = 2 * theta * oneMinusTheta; 97 final double coeffDot4 = -theta * oneMinus2Theta; 98 for (int i = 0; i < interpolatedState.length; ++i) { 99 final double yDot1 = yDotK[0][i]; 100 final double yDot23 = yDotK[1][i] + yDotK[2][i]; 101 final double yDot4 = yDotK[3][i]; 102 interpolatedState[i] = 103 currentState[i] + coeff1 * yDot1 + coeff23 * yDot23 + coeff4 * yDot4; 104 interpolatedDerivatives[i] = 105 coeffDot1 * yDot1 + coeffDot23 * yDot23 + coeffDot4 * yDot4; 106 } 107 108 } 109 110} 111