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.

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 -minDescendants-> 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.

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(t0) = y0$ the RK4 method is given by the following equations:

$$ \begin{align} k_1 &= y_n+\tfrac{1}{6}\left(k_1+2k_2+2k_3+k_4 \right)\\ t_{n+1} &= t_n+h \end{align} $$ where $y_{n + 1}$ is the RK4 approximation of $y(t_{n + 1})$, and

$$ \begin{align} k_1 &= f(t_n,y_n)\\ k_2 &= f(t_n+\tfrac{1}{2}h,y_n+\tfrac{1}{2}h k_1)\\ k_3 &= f(t_n+\tfrac{1}{2}h,y_n+\tfrac{1}{2}h k_2)\\ k_4 &= f(t_n+h,y_n+h k_3)\\ \end{align} $$

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

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 -minDescendants-> 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.

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.

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.

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

Trigger an event in regular intervals and call the event handler. The event handler must not modify the graph.

This allows to keep track of the changes for example:

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.addRow().setX(0, t).setY(0, x);
}
public void run ()
{
  monitorPeriodic(1, new Runnable() { 
    public void run() { data.addRow().setX(0, t).setY(0, x); }
  });

  integrate(10);
}
protected void getRate()
{
  t :'= 1;
  x :'= 0.1 * x;
}