core/org.eclipse.stem.solvers.rk/src/org/eclipse/stem/solvers/rk/impl/ODETest.java

package org.eclipse.stem.solvers.rk.impl;

import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.FirstOrderIntegrator;
import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
import org.apache.commons.math.ode.sampling.StepHandler;
import org.apache.commons.math.ode.sampling.StepInterpolator;

public class ODETest {

	public ODETest() {
		
		FirstOrderIntegrator dp853 = new DormandPrince853Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
		
		StepHandler stepHandler = new StepHandler() {
		    public void reset() {}
		            
		    public void handleStep(StepInterpolator interpolator, boolean isLast) throws DerivativeException {
		        double   t = interpolator.getCurrentTime();
		        double[] y = interpolator.getInterpolatedState();
		        System.out.println(t + " " + y[0] + " " + y[1]);
		    }

			public boolean requiresDenseOutput() {
				return true;
			}
		};
		dp853.addStepHandler(stepHandler);
		
		MyFunction ode = new ODETest.MyFunction(new double[] { 1.0, 1.0 }, 0.1);
		double[] y = new double[] { 0.0, 1.0 }; // initial state
		try {
			dp853.integrate(ode, 0.0, y, 16.0, y); // now y contains final state at time t=16.0
		} catch(Exception e) {
			e.printStackTrace();
		}
	}
	
	public static void main(String [] args) {
		
		new ODETest();
		
	}
	
	
	public class MyFunction implements FirstOrderDifferentialEquations {

		private double[] c;
	    private double omega;
	    
	    public MyFunction(double [] c, double omega) {
	    	this.c = c;
	    	this.omega = omega;
	    }
	    
		public int getDimension() {
			return 2;
		}

		public void computeDerivatives(double t, double[] y, double[] yDot)
				throws DerivativeException {
			 yDot[0] = omega * (c[1] - y[1]);
		       yDot[1] = omega * (y[0] - c[0]);
		}
		
	}
}