/*
* Open Source Physics software is free software as described near the bottom of this code file.
*
* For additional information and documentation on Open Source Physics please see:
* <http://www.opensourcephysics.org/>
*/
package org.opensourcephysics.numerics;
/**
* CashKarp45 implements a RKF 4/5 ODE solver with variable step size using Cash-Karp coefficients.
* @author W. Christian
* @author F. Esquembre
* @version 1.0
*/
public class CashKarp45 implements ODEAdaptiveSolver {
int error_code = ODEAdaptiveSolver.NO_ERROR;
// embedding constants Cash-Karp 4th and 5th order
static final double[][] a = {
{1.0/5.0}, {3.0/40.0, 9.0/40.0}, {3.0/10.0, -9.0/10.0, 6.0/5.0}, {-11.0/54.0, 5.0/2.0, -70.0/27.0, 35.0/27.0}, {1631.0/55296.0, 175.0/512.0, 575.0/13824.0, 44275.0/110592.0, 253.0/4096.0}
};
// b4 are 4th order coefficients
// static final double[] b4={ 2825./27648., 0., 18575./48384., 13525./55296., 277./14336., 1./4.};
// b5 are 5th oreer coefficients
static final double[] b5 = {37./378., 0., 250./621., 125./594., 0., 512./1771.};
static final double[] er = {277./64512., 0., -6925./370944., 6925./202752., 277./14336., -277./7084.};
static final int numStages = 6; // number of intermediate rate computations
private double stepSize = 0.01;
private int numEqn = 0;
private double[] temp_state;
private double[][] k;
private double truncErr;
private ODE ode;
protected double tol = 1.0e-6;
protected boolean enableExceptions = false;
/**
* Constructs the CashKarp45 ODESolver for a system of ordinary differential equations.
*
* @param _ode the system of differential equations.
*/
public CashKarp45(ODE _ode) {
ode = _ode;
initialize(stepSize);
}
/**
* Initializes the ODE solver.
*
* Temporary state and rate arrays are allocated.
* The number of differential equations is determined by invoking getState().length on the ODE.
*
* @param _stepSize
*/
public void initialize(double _stepSize) {
stepSize = _stepSize;
double state[] = ode.getState();
if(state==null) { // state vector not defined.
return;
}
if(numEqn!=state.length) {
numEqn = state.length;
temp_state = new double[numEqn];
k = new double[numStages][numEqn]; // six intermediate rates
}
}
/**
* Steps (advances) the differential equations by the stepSize.
*
* The ODESolver invokes the ODE's getRate method to obtain the initial state of the system.
* The ODESolver then advances the solution and copies the new state into the
* state array at the end of the solution step.
*
* @return the step size
*/
public double step() {
error_code = ODEAdaptiveSolver.NO_ERROR;
int iterations = 10;
double currentStep = stepSize, error = 0;
double state[] = ode.getState();
ode.getRate(state, k[0]); // get the initial rate
do {
iterations--;
currentStep = stepSize;
// Compute the k's
for(int s = 1; s<numStages; s++) {
for(int i = 0; i<numEqn; i++) {
temp_state[i] = state[i];
for(int j = 0; j<s; j++) {
temp_state[i] = temp_state[i]+stepSize*a[s-1][j]*k[j][i];
}
}
ode.getRate(temp_state, k[s]);
}
// Compute the error
error = 0;
for(int i = 0; i<numEqn; i++) {
truncErr = 0;
for(int s = 0; s<numStages; s++) {
truncErr = truncErr+stepSize*er[s]*k[s][i];
}
error = Math.max(error, Math.abs(truncErr));
}
if(error<=Float.MIN_VALUE) { // error too small to be meaningful,
error = tol/1.0e5; // increase stepSize x10
}
// find h step for the next try.
if(error>tol) { // shrink, no more than x10
double fac = 0.9*Math.pow(error/tol, -0.25);
stepSize = stepSize*Math.max(fac, 0.1);
} else if(error<tol/10.0) { // grow, but no more than factor of 10
double fac = 0.9*Math.pow(error/tol, -0.2);
if(fac>1) { // sometimes fac is <1 because error/tol is close to one
stepSize = stepSize*Math.min(fac, 10);
}
}
} while((error>tol)&&(iterations>0));
// advance the state
for(int i = 0; i<numEqn; i++) {
for(int s = 0; s<numStages; s++) {
state[i] += currentStep*b5[s]*k[s][i];
}
}
if(iterations==0) {
error_code = ODEAdaptiveSolver.DID_NOT_CONVERGE;
if(enableExceptions) {
throw new ODESolverException("CashKarp45 ODE solver did not converge."); //$NON-NLS-1$
}
}
return currentStep; // the value of the step actually taken.
}
/**
* Enables runtime exceptions if the solver does not converge.
* @param enable boolean
*/
public void enableRuntimeExpecptions(boolean enable) {
this.enableExceptions = enable;
}
/**
* Sets the step size.
*
* The step size may change when the step method is invoked.
*
* @param stepSize
*/
public void setStepSize(double stepSize) {
this.stepSize = stepSize;
}
/**
* Gets the step size.
*
* The stepsize is adaptive and may change as the step() method is invoked.
*
* @return the step size
*/
public double getStepSize() {
return stepSize;
}
/**
* Method setTolerance
*
* @param _tol
*/
public void setTolerance(double _tol) {
tol = Math.abs(_tol);
if(tol<1.0E-12) {
String err_msg = "Error: CashKarp ODE solver tolerance cannot be smaller than 1.0e-12."; //$NON-NLS-1$
if(enableExceptions) {
throw new ODESolverException(err_msg);
}
System.err.println(err_msg);
}
}
/**
* Method getTolerance
*
* @return tolerance
*/
public double getTolerance() {
return tol;
}
/**
* Gets the error code.
* Error codes:
* ODEAdaptiveSolver.NO_ERROR
* ODEAdaptiveSolver.DID_NOT_CONVERGE
* @return int
*/
public int getErrorCode() {
return error_code;
}
}
/*
* Open Source Physics software is free software; you can redistribute
* it and/or modify it under the terms of the GNU General Public License (GPL) as
* published by the Free Software Foundation; either version 2 of the License,
* or(at your option) any later version.
* Code that uses any portion of the code in the org.opensourcephysics package
* or any subpackage (subdirectory) of this package must must also be be released
* under the GNU GPL license.
*
* This software is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
* or view the license online at http://www.gnu.org/copyleft/gpl.html
*
* Copyright (c) 2007 The Open Source Physics project
* http://www.opensourcephysics.org
*/