From groimp

Since GroIMP 1.0 a new ODE framework is available to modellers. This allows easy specification and solution of differential equations, which are used in many FSPMs. The key component of this framework is a new operator :'= called the deferred rate assignment operator.

Contents 1 Introductory example 1.1 Advanced numerical methods 1.2 Using the ODE framework 1.3 The getRate() function 2 Simple example 3 Monitor functions 4 API 5 FAQ

Introductory example

Let's start with a simple example. To simulate transport of a substrate (carbon assimilates) through a structure (connected organs of a virtual plant), diffusion between neighbouring parts of the plant is modelled. Without the ODE framework, this could (and in fact was, by some users) written as XL code as follows:

a:C -neighbor-> b:C ::> {
  double rate = D * (b[carbon] - a[carbon]);
  a[carbon] :+= h*rate;
  b[carbon] :-= h*rate;
}

From a numerical point of view this corresponds to simple explicit Euler integration (with step size of h, often just set to 1). This should be avoided in practice, since the results from an Euler numerical integration can differ substantially from the exact analytical solution, depending on the type of ODE and the integration step size chosen. Euler integration also suffers from stability problems if the step size chosen is too large.

Advanced numerical methods

A very commonly used method to solve initial value problems is the classical Runge-Kutta method of 4th order. This method is simple, fast, and robust.

For an initial value problem y' = f(t,y),y(t₀) = y₀ the RK4 method is given by the following equations:

where y_{n + 1} is the RK4 approximation of y(t_{n + 1}), and

are the values of the first derivation for different combinations of time and state.

Using the ODE framework

The interface between the problem to be solved and an ODE integrator like RK4 is given by the initial state and time and the rate function f(t,y). Specification of ODEs from within an RGG must occur inside a function void getRate(). The above example of diffusion then would be rewritten to use the operator :'= as follows:

a:C -neighbor-> b:C ::> {
  double rate = D * (b[carbon] - a[carbon]);
  a[carbon] :'= +rate;
  b[carbon] :'= -rate;
}

Note that the integration step size is not given explicitly, as choosing it is up to the integrator. While the two ways of defining the diffusion of the carbon assimilates look very similar (which was an important design goal), both of them are handled very differently internally.

The getRate() function

All ODEs must be specified by the user inside a function getRate(). The current state is provided implicitly as the attributes of the graph. The user should not perform structural changes to the graph like adding/removing nodes. This will result in undefined behavior. If structural changes need to be done, use a monitor that stops integration and perform the structural change there.

Simple example

A simple but complete example as RGG program is this:

module A(double len);
protected void init ()
[
  Axiom ==> A(1);
]
public void run ()
{
  integrate(10);
  println((* a:A *)[len] + " (" + a + ")");
}
protected void getRate()
[
  a:A ::> a[len] :'= 1.23;
]

Here, a module A is created with an initial value of its attribute len of 1. Integration is performed over 10 time units. Afterwards the resulting value of the attribute is printed. The rate function f(t,y) is given by the getRate() function and defines a slope of 1.23 for the len attribute.

Monitor functions

A common situation that might appear is to stop an integration conditionally, i.e. once the output of the integration has reached a certain value. Consider the simple example from above. If we would like to know the time when len has reached a value of 10 we somehow have to interact with the integrator. The concept used here is called monitor function, sometimes also called trigger function. The user provides a function g(y) that calculates a single value for a given state. During integration this function is evaluated to detect sign changes. When this happens, root finding algorithms can be used to find the exact time when the sign change occurred. The modified example then looks like this:

module A(double len);
protected void init ()
[
  Axiom ==> A(1);
]
public void run ()
{
  [ a:A ::> monitor(void=>double a[len] - 10); ]
  integrate(10);
  println((* a:A *)[len] + " (" + a + ")");
}
protected void getRate()
[
  a:A ::> a[len] :'= 1.23;
]

The monitor() function is provided as part of the public API of the ODE framework and expects a functor as its parameter representing the function g(y), where the state y is provided implicitly as the state of the graph when the functor is executed. Here we just calculate a[len] - 10, so that when the attribute reaches a value of 10 integration will be stopped.

Together with the monitor function, also an event handler may be specified. This will be executed once the correct time for the event was found. Imagine you want to plot the state at regular intervals. An example program is this:

import static java.lang.Math.*;
DatasetRef data = new DatasetRef("Dataplot");
const double PERIOD = 1.0;
double t;
double x;
protected void init ()
{
  data.clear().setColumnKey(0, "x");
  chart(data, XY_PLOT);
  t = 0;
  x = 1;
  data << x;
}
public void run ()
{
  monitor(void=>double sin(PI*t/PERIOD), new Runnable() { 
    public void run() { data << x; }
  }, CONTINUE, 0.25*PERIOD, 1e-4, 100);
  integrate(20);
}
protected void getRate()
{
  t :'= 1;
  x :'= 0.1 * x;
}

Note that the values t and x as variables of RGG can be integrated as well, because the RGG itself is also a node. In general only attributes of nodes may be integrated. Besides the monitor function and the event handle, some additional parameters are provided in the call to monitor(). The flag CONTINUE tells the integrator to execute the event and then continue integration irrespective of previous events. This is available for cases like those shown above, where the event handler wants to extract the current state (or parts of it). The state of the graph should not be modified when using this flag. The value of 0.25*PERIOD is used to tell the integrator with which maximal distance the monitor function should be evaluated. This duration must be chosen small enough so that no sign change is missed, but should be chosen as big as possible to reduce computational load. The value of 1e-4 defines the accuracy of the root to look for and the value of 100 gives the number of iterations.

API

void integrate(double duration)

Perform integration of the ODEs for the given duration. When monitor functions are set, the function may also return before duration is reached.

void integrate()

This just calls integrate(Double.MAX_VALUE). A monitor function must have been set to stop integration.

void monitor(final VoidToDouble g)
void monitor(final VoidToDouble g, final Runnable r)
void monitor(final VoidToDouble g, final VoidToVoid vtv)

Set a monitor that stops integration when the sign changes. This can be used, for instance, to integrate until a value reaches a threshold. The monitor function is a functor defined by preceding an expression by void=>double. The result type of that expression must be of type double. If a second parameter is given, this is used as event handler and is executed at the moment g changes its sign.

void monitor(final VoidToDouble g, final Runnable r, final int action, final double maxCheckInterval,
           final double convergence, final int maxIterationCount)

Set a monitor that triggers when g changes its sign. Root finding is used to obtain the time when the sign change occurs within an accuracy of convergence time units and a maximum iteration count for the root finding algorithm of maxIterationCount. After the event handler r has been executed, action determines what to do next:

• STOP stop integration at this point
• CONTINUE continue integration, no matter what happened before
• RESET_STATE informs the integrator that the current state has changed and thus is discontinuous
• RESET_DERIVATIVES informs the integrator that the derivative is discontinuous

(Note: the RESET_* actions are not supported yet)

Solver getSolver()

Returns the currently set solver. By default, a DormandPrince54Integrator is set. Its minimum step is zero, its maximum step is one, and its relative and absolute error is 1e-4.

void setSolver(Solver solver)
void setSolver(FirstOrderIntegrator solver)

Set a new solver as current solver. In case the default solver causes problems, for instance, because special parameters must be chosen, this allows to choose a more appropriate solver. More information about possible solvers can be found here. For instance you can call setSolver(new org.apache.commons.math.ode.nonstiff.GraggBulirschStoerIntegrator(0, 1, 1E-4, 1E-4)) to change the current solver to GBS with minimum step size of zero, maximum step size of one, and a maximum allowed relative and absolute error of 1E-4.

FAQ

• Is there an upper limit to the number of rates that can be used in the getRate() method?
  

Besides the obvious limit of available system memory, there can be at most 2³¹ − 1 values being integrated. This is because internally the integrator uses double arrays to store state and rate vectors and Java permits only a signed 32-bit integer value as size when allocating the array. Other restrictions may apply depending on the internal implementation of the used integrator.

• The rates that you show in your examples look rather simple. I am using a photosynthesis model which outputs instantaneous rates of net assimilation as a function of several environment variables. This is an external Java file. Can the ODE solver handle this?
  

Yes, just call the photosynthesis model from within your getRate() function and apply the rates using :'=.

• In the photosynthesis model mentioned above the environmental input data are at different, yet fixed resolutions (10 minutes, 1 hour). With the Euler integration I have used before this was not a problem but you mention the advanced RK solvers choose their own step size.
  

If you have measured input data you can interpolate them to make them continuous, then feed that continuous data to the photosynthesis model to calculate rates. Interpolation can be done using SplineFunction, for instance.
Regarding the use of variable step sizes, this depends on the implementation of the integrator. The integrators provided by Apache Commons Math are listed here. For instance, the default integrator Dormand-Prince 5(4) provides variable step size control, while the Classical Runge-Kutta uses a fixed-size step.

• What is Apache Commons Math ?
  

The Commons is an Apache project focused on all aspects of reusable Java components. (from their website)
Apache Commons Math is part of Apache Commons and is a library dedicated to mathematics and statistics. A part of that math library are numerical integrators, that can be used as described above.

• I set ClassicalRungeKuttaIntegrator as solver, but when I call getSolver() the system returns me de.grogra.numeric.FirstOrderIntegratorAdapter. Is that correct ?
  

Yes, this is correct. All solvers provided by an external library like Apache Commons Math will have to be wrapped inside another class, that performs the necessary work to make the compatible to the ODE system. In this case the wrapper class FirstOrderIntegratorAdapter allows to wrap solvers that implement FirstOrderIntegrator, such as those provided by Apache Commons Math.

write your questions here