ODEAdaptiveSolver

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;
  }
}