ODESolverRK4

package odend;

import Jama.*;

/** 
 * solver using classical 4th order Runge-Kutta method
 */
public class ODESolverRK4 extends ODESolver {

  /** sets ode and initial values */
  public ODESolverRK4(ODE ode) {
    super(ode);
  }

  /** 
   * computes solution at next time step<br>
   * using the classical 4th order Runke-Kutta scheme
   */
  public void nextStep(double h) {
    Matrix k1 = ode.f(x, t);
    
    double t1 = t + h / 2.0;
    Matrix x1 = x.plus(k1.times(h / 2.0)); // x1 = x + h*k1/2.0
    Matrix k2 = ode.f(x1, t1);
    
    double t2 = t1;	
    Matrix x2 = x.plus(k2.times(h / 2.0)); // x2 = x + h*k2/2.0
    Matrix k3 = ode.f(x2, t2);
    
    double t3 = t + h;
    Matrix x3 = x.plus(k3.times(h)); // x3 = x + k3*h
    Matrix k4 = ode.f(x3, t3);
    
    // x += (k1 + 2.0*k2 + 2.0*k3 + k4)*h/6.0
    x.plusEquals(k1.times(h / 6.0));
    x.plusEquals(k2.times(h / 3.0));
    x.plusEquals(k3.times(h / 3.0));
    x.plusEquals(k4.times(h / 6.0));
    
    t += h;
  }
}