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