package odeadapt; import Jama.*; /** * adaptive solver that uses a special strategy to get a new step size<br> * <br> * basic strategy: <br> * integrate with actual stepsize in one step<br> * integrate with half of actual stepsize in two steps<br> * estimate the error to get new stepsize */ public class ODEAdaptiveSolver extends ODESolver { protected double DEFAULT_TOLERANCE = 1.0/1000000.0; protected ODESingleStepSolver solv; protected int order; // order of solver = order of ode1 // end point of integration protected double tEnd; // current value of step size protected double stepSize; // max. error protected double tol; // factor for error estimation protected double gamma; // safety factor for step size protected double safety = 95.0/100.0; // maximal expansion factor of step size protected double maxFac = 5.0; /** * construct adaptive solver for a given ODE defining<br> * end time of integration<br> * maximal local error<br> * a single step solver for the integration in one and two half steps */ public ODEAdaptiveSolver(ODE ode, double tEnd, ODESingleStepSolver solv) { super(ode); // sets t, x this.order = solv.order; this.solv = solv; if (tEnd > ode.t0) { this.tEnd = tEnd; } else { // tEnd is smaller than t0! throw new IllegalArgumentException("Start time " + ode.t0 + " is later than end time " + tEnd); } // try to integrate in one step at first stepSize = tEnd - ode.t0; gamma = 1.0 / (1.0 - Math.pow(2.0, -order)); } /** * construct adaptive solver for a given ODE defining<br> * end time of integration<br> * maximal local error<br> * use Runge-Kutta-4 solver as single step solver */ public ODEAdaptiveSolver(ODE ode, double tEnd) { this(ode, tEnd, new ODESolverRK4(ode)); } /** * integrate until t+dt<br> * return number of steps needed */ public int nextStep(double dt) throws StepsizeTooSmallException { // old stepSize could be too small while approaching tEnd from last step! // therefore: try with one step // alternatively: use standard initial step size stepSize = dt; int nSteps = 0; tEnd = t + dt; while (t < tEnd - tol) { nextStep(); nSteps++; } return nSteps; } /** * proceed one step, using an appropriate step size */ public void nextStep() throws StepsizeTooSmallException { boolean stepAccepted = false; while (!stepAccepted) { // solve in one large step solv.t = t; solv.x = x.copy(); solv.nextStep(stepSize); Matrix xOneStep = solv.x; // solve in two small steps solv.t = t; solv.x = x.copy(); solv.nextStep(stepSize / 2.0); solv.nextStep(stepSize / 2.0); // estimate error // diff = gamma * (x2 - x1) Matrix diff = solv.x.minus(xOneStep).timesEquals(gamma); double error = diff.normF(); // compute new stepsize double hRel = Math.pow(tol / error, 1.0 / (1.0 + order)); double stepSizeNew = stepSize * hRel; // throw exception, if step size gets too small if (stepSizeNew < tol) { throw new StepsizeTooSmallException("Step size " + stepSize + " is smaller than tolerance " + tol); } if (hRel < 1) { // step rejected, use smaller step size stepSize = safety * stepSizeNew; } else { // step accepted, use new values stepAccepted = true; t += stepSize; x = xOneStep.plus(diff); // increase step size if (hRel > maxFac) { // don't expand too much ! stepSize *= maxFac; } else { stepSize = safety * stepSizeNew; } stepSize = Math.min(stepSize, tEnd - t); // don't overshoot! } } } /** * set absolute tolerance<br> * default usually: 1e-6 */ public void setAbsoluteTolerance(double t) { tol = t; } /** * set maximal factor m<br> * new stepsize is smaller than m *old_stepsize <br> * should be larger than 1 */ public void setMaxFac(double m) { maxFac = m; } /** * set safety factor s<br> * multiplies proposed new stepsize<br> * should be smaller than 1 <br> * usually approx. 1 */ public void setSafety(double s) { safety = s; } }