Explicit fixed-step methods

Functions
array[] vector	odeint_euler (real t0, vector y0, data real h, data int num_steps, array[] int a0, vector theta)
array[] vector	odeint_midpoint (real t0, vector y0, data real h, data int num_steps, array[] int a0, vector theta)
array[] vector	odeint_rk4 (real t0, vector y0, data real h, data int num_steps, array[] int a0, vector theta)

Detailed Description

These are explicit fixed-step ODE integrators, which evaluate the solution at an equispaced grid of time points. Interpolation can then be used to get the solution at other desired output time points. The functions assume that derivative_fun(), which is the system function of the ODE, is defined earlier in the functions block. It must have signature

vector derivative_fun(real t, vector y, int[] a0, vector a1)

i.e. same as what can be used with the built-in ODE solvers of Stan.

 
array[] vector odeint_euler(real t0, vector y0, data real h,
                            data int num_steps, array[] int a0, vector theta) {
  int d = num_elements(y0);
  array[num_steps + 1] vector[d] y;
  real t = t0;
  y[1] = y0;
  
  // Integrate at grid of time points
  for (i in 1 : num_steps) {
    y[i + 1] = y[i] + h * derivative_fun(t, y[i], a0, theta);
    t = t + h;
  }
  return y;
}
 

 
array[] vector odeint_midpoint(real t0, vector y0, data real h,
                               data int num_steps, array[] int a0,
                               vector theta) {
  int d = num_elements(y0);
  array[num_steps + 1] vector[d] y;
  vector[d] y_mid;
  real t = t0;
  y[1] = y0;
  
  // Integrate at grid of time points
  for (i in 1 : num_steps) {
    // Half-Euler step
    y_mid = y[i] + 0.5 * h * derivative_fun(t, y[i], a0, theta);
    
    // Full step using derivative at midpoint
    y[i + 1] = y[i] + h * derivative_fun(t + 0.5 * h, y_mid, a0, theta);
    t = t + h;
  }
  
  return y;
}

 
array[] vector odeint_rk4(real t0, vector y0, data real h,
                          data int num_steps, array[] int a0, vector theta) {
  int d = num_elements(y0);
  array[num_steps + 1] vector[d] y;
  vector[d] k1;
  vector[d] k2;
  vector[d] k3;
  vector[d] k4;
  real t = t0;
  y[1] = y0;
  
  // Integrate at grid of time points
  for (i in 1 : num_steps) {
    k1 = h * derivative_fun(t, y[i], a0, theta);
    k2 = h * derivative_fun(t + 0.5 * h, y[i] + 0.5 * k1, a0, theta);
    k3 = h * derivative_fun(t + 0.5 * h, y[i] + 0.5 * k2, a0, theta);
    k4 = h * derivative_fun(t + h, y[i] + k3, a0, theta);
    y[i + 1] = y[i] + (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;
    t = t + h;
  }
  
  return y;
}

Function Documentation

odeint_euler()

array[] vector odeint_euler	(	real	t0,
		vector	y0,
		data real	h,
		data int	num_steps,
		array[]int	a0,
		vector	theta
	)

Forward Euler method

Info: https://en.wikipedia.org/wiki/Euler_method

Author

Juho Timonen

Parameters

t0	initial time
y0	initial state, D-vector
h	step size (positive)
num_steps	number of steps to take
a0	array of integer inputs given to derivative_fun()
theta	parameter vector given to derivative_fun()

Returns

array of D-vectors, length equal to num_steps + 1

odeint_midpoint()

array[] vector odeint_midpoint	(	real	t0,
		vector	y0,
		data real	h,
		data int	num_steps,
		array[]int	a0,
		vector	theta
	)

Explicit midpoint method

Info: https://en.wikipedia.org/wiki/Midpoint_method

Author

Juho Timonen

Parameters

t0	initial time
y0	initial state, D-vector
h	step size (positive)
num_steps	number of steps to take
a0	array of integer inputs given to derivative_fun()
theta	parameter vector given to derivative_fun()

Returns

array of D-vectors, length equal to num_steps + 1

odeint_rk4()

array[] vector odeint_rk4	(	real	t0,
		vector	y0,
		data real	h,
		data int	num_steps,
		array[]int	a0,
		vector	theta
	)

Fourth-order Runge-Kutta method

Info: https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

Author

Juho Timonen

Parameters

t0	initial time
y0	initial state, D-vector
h	step size (positive)
num_steps	number of steps to take
a0	array of integer inputs given to derivative_fun()
theta	parameter vector given to derivative_fun()

Returns

array of D-vectors, length equal to num_steps + 1