001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018package org.apache.commons.math4.legacy.ode.sampling; 019 020import java.io.IOException; 021import java.io.ObjectInput; 022import java.io.ObjectOutput; 023import java.util.Arrays; 024 025import org.apache.commons.math4.legacy.exception.MaxCountExceededException; 026import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix; 027import org.apache.commons.math4.legacy.ode.EquationsMapper; 028import org.apache.commons.math4.core.jdkmath.JdkMath; 029 030/** 031 * This class implements an interpolator for integrators using Nordsieck representation. 032 * 033 * <p>This interpolator computes dense output around the current point. 034 * The interpolation equation is based on Taylor series formulas. 035 * 036 * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsBashforthIntegrator 037 * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsMoultonIntegrator 038 * @since 2.0 039 */ 040 041public class NordsieckStepInterpolator extends AbstractStepInterpolator { 042 043 /** Serializable version identifier. */ 044 private static final long serialVersionUID = -7179861704951334960L; 045 046 /** State variation. */ 047 protected double[] stateVariation; 048 049 /** Step size used in the first scaled derivative and Nordsieck vector. */ 050 private double scalingH; 051 052 /** Reference time for all arrays. 053 * <p>Sometimes, the reference time is the same as previousTime, 054 * sometimes it is the same as currentTime, so we use a separate 055 * field to avoid any confusion. 056 * </p> 057 */ 058 private double referenceTime; 059 060 /** First scaled derivative. */ 061 private double[] scaled; 062 063 /** Nordsieck vector. */ 064 private Array2DRowRealMatrix nordsieck; 065 066 /** Simple constructor. 067 * This constructor builds an instance that is not usable yet, the 068 * {@link AbstractStepInterpolator#reinitialize} method should be called 069 * before using the instance in order to initialize the internal arrays. This 070 * constructor is used only in order to delay the initialization in 071 * some cases. 072 */ 073 public NordsieckStepInterpolator() { 074 } 075 076 /** Copy constructor. 077 * @param interpolator interpolator to copy from. The copy is a deep 078 * copy: its arrays are separated from the original arrays of the 079 * instance 080 */ 081 public NordsieckStepInterpolator(final NordsieckStepInterpolator interpolator) { 082 super(interpolator); 083 scalingH = interpolator.scalingH; 084 referenceTime = interpolator.referenceTime; 085 if (interpolator.scaled != null) { 086 scaled = interpolator.scaled.clone(); 087 } 088 if (interpolator.nordsieck != null) { 089 nordsieck = new Array2DRowRealMatrix(interpolator.nordsieck.getDataRef(), true); 090 } 091 if (interpolator.stateVariation != null) { 092 stateVariation = interpolator.stateVariation.clone(); 093 } 094 } 095 096 /** {@inheritDoc} */ 097 @Override 098 protected StepInterpolator doCopy() { 099 return new NordsieckStepInterpolator(this); 100 } 101 102 /** Reinitialize the instance. 103 * <p>Beware that all arrays <em>must</em> be references to integrator 104 * arrays, in order to ensure proper update without copy.</p> 105 * @param y reference to the integrator array holding the state at 106 * the end of the step 107 * @param forward integration direction indicator 108 * @param primaryMapper equations mapper for the primary equations set 109 * @param secondaryMappers equations mappers for the secondary equations sets 110 */ 111 @Override 112 public void reinitialize(final double[] y, final boolean forward, 113 final EquationsMapper primaryMapper, 114 final EquationsMapper[] secondaryMappers) { 115 super.reinitialize(y, forward, primaryMapper, secondaryMappers); 116 stateVariation = new double[y.length]; 117 } 118 119 /** Reinitialize the instance. 120 * <p>Beware that all arrays <em>must</em> be references to integrator 121 * arrays, in order to ensure proper update without copy.</p> 122 * @param time time at which all arrays are defined 123 * @param stepSize step size used in the scaled and Nordsieck arrays 124 * @param scaledDerivative reference to the integrator array holding the first 125 * scaled derivative 126 * @param nordsieckVector reference to the integrator matrix holding the 127 * Nordsieck vector 128 */ 129 public void reinitialize(final double time, final double stepSize, 130 final double[] scaledDerivative, 131 final Array2DRowRealMatrix nordsieckVector) { 132 this.referenceTime = time; 133 this.scalingH = stepSize; 134 this.scaled = scaledDerivative; 135 this.nordsieck = nordsieckVector; 136 137 // make sure the state and derivatives will depend on the new arrays 138 setInterpolatedTime(getInterpolatedTime()); 139 } 140 141 /** Rescale the instance. 142 * <p>Since the scaled and Nordsieck arrays are shared with the caller, 143 * this method has the side effect of rescaling this arrays in the caller too.</p> 144 * @param stepSize new step size to use in the scaled and Nordsieck arrays 145 */ 146 public void rescale(final double stepSize) { 147 148 final double ratio = stepSize / scalingH; 149 for (int i = 0; i < scaled.length; ++i) { 150 scaled[i] *= ratio; 151 } 152 153 final double[][] nData = nordsieck.getDataRef(); 154 double power = ratio; 155 for (int i = 0; i < nData.length; ++i) { 156 power *= ratio; 157 final double[] nDataI = nData[i]; 158 for (int j = 0; j < nDataI.length; ++j) { 159 nDataI[j] *= power; 160 } 161 } 162 163 scalingH = stepSize; 164 } 165 166 /** 167 * Get the state vector variation from current to interpolated state. 168 * <p>This method is aimed at computing y(t<sub>interpolation</sub>) 169 * -y(t<sub>current</sub>) accurately by avoiding the cancellation errors 170 * that would occur if the subtraction were performed explicitly.</p> 171 * <p>The returned vector is a reference to a reused array, so 172 * it should not be modified and it should be copied if it needs 173 * to be preserved across several calls.</p> 174 * @return state vector at time {@link #getInterpolatedTime} 175 * @see #getInterpolatedDerivatives() 176 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 177 */ 178 public double[] getInterpolatedStateVariation() throws MaxCountExceededException { 179 // compute and ignore interpolated state 180 // to make sure state variation is computed as a side effect 181 getInterpolatedState(); 182 return stateVariation; 183 } 184 185 /** {@inheritDoc} */ 186 @Override 187 protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) { 188 189 final double x = interpolatedTime - referenceTime; 190 final double normalizedAbscissa = x / scalingH; 191 192 Arrays.fill(stateVariation, 0.0); 193 Arrays.fill(interpolatedDerivatives, 0.0); 194 195 // apply Taylor formula from high order to low order, 196 // for the sake of numerical accuracy 197 final double[][] nData = nordsieck.getDataRef(); 198 for (int i = nData.length - 1; i >= 0; --i) { 199 final int order = i + 2; 200 final double[] nDataI = nData[i]; 201 final double power = JdkMath.pow(normalizedAbscissa, order); 202 for (int j = 0; j < nDataI.length; ++j) { 203 final double d = nDataI[j] * power; 204 stateVariation[j] += d; 205 interpolatedDerivatives[j] += order * d; 206 } 207 } 208 209 for (int j = 0; j < currentState.length; ++j) { 210 stateVariation[j] += scaled[j] * normalizedAbscissa; 211 interpolatedState[j] = currentState[j] + stateVariation[j]; 212 interpolatedDerivatives[j] = 213 (interpolatedDerivatives[j] + scaled[j] * normalizedAbscissa) / x; 214 } 215 } 216 217 /** {@inheritDoc} */ 218 @Override 219 public void writeExternal(final ObjectOutput out) 220 throws IOException { 221 222 // save the state of the base class 223 writeBaseExternal(out); 224 225 // save the local attributes 226 out.writeDouble(scalingH); 227 out.writeDouble(referenceTime); 228 229 final int n = (currentState == null) ? -1 : currentState.length; 230 if (scaled == null) { 231 out.writeBoolean(false); 232 } else { 233 out.writeBoolean(true); 234 for (int j = 0; j < n; ++j) { 235 out.writeDouble(scaled[j]); 236 } 237 } 238 239 if (nordsieck == null) { 240 out.writeBoolean(false); 241 } else { 242 out.writeBoolean(true); 243 out.writeObject(nordsieck); 244 } 245 246 // we don't save state variation, it will be recomputed 247 } 248 249 /** {@inheritDoc} */ 250 @Override 251 public void readExternal(final ObjectInput in) 252 throws IOException, ClassNotFoundException { 253 254 // read the base class 255 final double t = readBaseExternal(in); 256 257 // read the local attributes 258 scalingH = in.readDouble(); 259 referenceTime = in.readDouble(); 260 261 final int n = (currentState == null) ? -1 : currentState.length; 262 final boolean hasScaled = in.readBoolean(); 263 if (hasScaled) { 264 scaled = new double[n]; 265 for (int j = 0; j < n; ++j) { 266 scaled[j] = in.readDouble(); 267 } 268 } else { 269 scaled = null; 270 } 271 272 final boolean hasNordsieck = in.readBoolean(); 273 if (hasNordsieck) { 274 nordsieck = (Array2DRowRealMatrix) in.readObject(); 275 } else { 276 nordsieck = null; 277 } 278 279 if (hasScaled && hasNordsieck) { 280 // we can now set the interpolated time and state 281 stateVariation = new double[n]; 282 setInterpolatedTime(t); 283 } else { 284 stateVariation = null; 285 } 286 } 287}