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