Wednesday 23 July 2014

ODEs with Java using the Apache Commons Math library

It drives me nuts when I struggle to get my head around something just to, when I eventually do, discover that it is not really as complicated as the documentation makes it out to be. Simple examples are what I need. Here is an example of how to implement the Lorenz System. Not that I really know much about the Lorenz System but it is always a simple system to use when you want to try out a new way of solving differential equations. More information about the Lorenz System can be found in Wikipedia.

The first thing to do is to create a class containing the equations. The Lorenz system consists of the following three equations:


Doing this using the Apache Commons Math library means that we have to create a class that implements the FirstOrderDifferentialEquations interface. This interface requires two methods to be implemented, which are computeDerivatives and getDimension:


01 import org.apache.commons.math3.exception.DimensionMismatchException;
02 import org.apache.commons.math3.exception.MaxCountExceededException;
03 import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
04 
05 public class Lorenz implements FirstOrderDifferentialEquations {
06 
07         double sigma = 10.0;
08         double beta = 3;
09         double rho = 28.0;
10         double X, Y, Z;
11 
12         @Override
13         public void computeDerivatives(double t, double[] y, double[] yDot)
14                         throws MaxCountExceededException, 

                           DimensionMismatchException {
15                 X = y[0];
16                 Y = y[1];
17                 Z = y[2];
18                 yDot[0= sigma * (Y - X);
19                 yDot[1= X * (rho - Z- Y;
20                 yDot[2(X * Y(beta * Z);
21         }
22 
23         @Override
24         public int getDimension() {
25                 // TODO Auto-generated method stub
26                 return 3;
27         }
28 }
Java2html

The next step is to write the code that will call the integrator. There are various integrators. For this application we'll use the ClassicalRungeKuttaIntegrator. You will also need a step handler to capture each step of the integration:

01 import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
02 import org.apache.commons.math3.ode.nonstiff.ClassicalRungeKuttaIntegrator;
03 import org.apache.commons.math3.ode.sampling.StepHandler;
04 
05 public class LorenzMain {
06 
07         public static void main(String[] args) {
08 
09                 ClassicalRungeKuttaIntegrator rk4 = 

                           new ClassicalRungeKuttaIntegrator(0.01);
10                 FirstOrderDifferentialEquations ode = new Lorenz();
11                 double[] y = new double[] { 10.0, -2.050 }// initial state
12                 StepHandler stepHandler = new 

                           WriteToFileStepHandler("lorenz.csv");
13                 rk4.addStepHandler(stepHandler);
14                 rk4.integrate(ode, // equations
15                                 0.0// start time
16                                 y, // initial conditions
17                                 100.0// end time
18                                 y)// result
19                 for (int i = 0; i < y.length; i++) {
20                         System.out.println(y[i]);
21                 }
22         }
23 }
Java2html

The step handler has to implement the StepHandler class. For the step handler below I created a constructor that takes a filename so that the results can be written to a csv file. Obviously you can do whatever you want with this data:

01 import java.io.File;
02 import java.io.PrintWriter;
03 import java.util.ArrayList;
04 
05 import org.apache.commons.math3.exception.MaxCountExceededException;
06 import org.apache.commons.math3.ode.sampling.StepHandler;
07 import org.apache.commons.math3.ode.sampling.StepInterpolator;
08 
09 
10 public class WriteToFileStepHandler implements StepHandler {
11         ArrayList<String> steps = new ArrayList<String>();
12         String filename;
13         
14         public WriteToFileStepHandler(String filename) {
15                 this.filename = filename;
16         }
17         @Override
18         public void handleStep(StepInterpolator interpolator, boolean isLast)
19                         throws MaxCountExceededException {
20                 double t = interpolator.getCurrentTime();
21                 double[] y = interpolator.getInterpolatedState();
22                         String onestep = "" + t;
23                         for (int i = 0; i < y.length; i++) {
24                                 onestep += ", " + y[i];
25                         }
26                         steps.add(onestep);
27                 if (isLast) {
28                         try {
29                                 PrintWriter writer = 

30                                 new PrintWriter(new File(filename)
                                   "UTF-8");                                    
31                                 for (String step : steps) {
32                                         writer.println(step);
33                                 }
34                                 writer.close();
35                         catch (Exception e) {
36                         }
37                         ;
38                 }
39         }
40 
41         @Override
42         public void init(double t0, double[] y0, double t) {
43                 // TODO Auto-generated method stub
44 
45         }
46 
47 }
Java2html
Last but not the least I read the file with R and created the graph: