Re: Calling lsode directly from a c++ program via liboctave?

[Douglas Eck]

I figured it out. I think it's kind of cool. And fast.
Attached is an example that shows an example
of calling the ODE solver lsode directly in
a dynamically-linked c++ file as opposed to calling
it from octave.

The example used is a stiff ODE in the form of a
Fitzhugh-Nagumo relaxation oscillator. This oscillator
is a model of a neuron being driven with constant voltage
from a "space clamp" on it's axon. It is a two-varible
simplification of the famous Hodgkin Huxley equations.

If anyone can improve on it, please let me know.

Cheers,
Doug

fitz_driver.m
Description: type

/*
To make, type mkoctfile fitz.cc

This was tested using octave 2.1.43. I have
no idea if it works in 2.0

The example used is a stiff ODE in the form of a 
Fitzhugh-Nagumo relaxation oscillator. This oscillator
is a model of a neuron being driven with constant voltage
from a "space clamp" on it's axon. It is a two-varible
simplification of the famous Hodgkin Huxley equations. 

The Fitzhugh-Nagumo equations were derived simultaneoulsly by Fitzhugh and 
Nagumo et.al.

R. Fitzhugh, Impulses and physiological states in theoretical 
models of nerve membranes, 1182-Biophys. J., 1 (1961), pp. 445--466. and

J. Nagumo, S. Arimoto, and S. Yoshizawa, An active pulse 
transmission line simulating 1214-nerve axons, Proc. IRL, 50 (1960), pp. 
2061--2070.

Douglas Eck address@hidden

This code borrows heavily from John Eaton's lsode.cc
*/

#include <octave/config.h>
#include <octave/LSODE.h>
#include <octave/variables.h>
#include <octave/oct.h>
#include <octave/ov-struct.h>

#define FITZ_THOLD .2
#define FITZ_CVOLT .112
#define FITZ_SHUNT 1.2

static LSODE_options lsode_opts;
static ColumnVector fn_ode_epsilons;

ColumnVector fn_ode_helper (const ColumnVector& x, double t)
{
  //This is the function called multiple times by lsode.
  //Each oscillator has a voltage v and recovery w
  //thus x is of length i=oscCount*2 with x(i)=v, x(i+1)=w
  //The parameter epsilon controls the rate of oscillaton
  //it is stored in a static column vector fn_ode_epsilons
  int l=x.length();
  ColumnVector xdot(l);
  xdot.fill(0.0);
  for (int i=0;i<l;i=i+2) {
    xdot(i) = -x(i)*(x(i)-FITZ_THOLD) * (x(i)-1) - x(i+1) + FITZ_CVOLT;
    xdot(i+1) = fn_ode_epsilons(i/2) * (x(i)-FITZ_SHUNT * x(i+1));
  }
  return xdot;
}

DEFUN_DLD (fitz, args,,
"Fitzhugh-Nagumo oscillator solved using lsode. 

Calling syntax: fitz(x0,epsilons,t) where x0
is a vector defining starting conditions 
of the oscillators. Each oscillator is represented
using a duple (v,w)  v=starting voltage
 w=starting recovery, with initial conditions in
x0=(v1,w1,v2,w2,...,vK,wK) where K=oscCount

Epsilon is the variable controlling the rate (period) of the
oscillator. It is stored in a static ColumnVector.
")
{
  octave_value_list retval;
  if (args.length()>2 | args.length()<2) {
    print_usage("fitz");
  }
  ColumnVector state (args(0).vector_value ());
  fn_ode_epsilons=ColumnVector(args(1).vector_value());
  ColumnVector out_times (args(2).vector_value ());
  double tzero = out_times (0);
  int nsteps = out_times.capacity ();
  ODEFunc func (fn_ode_helper);
  LSODE ode (state, tzero, func);  
  ode.copy (lsode_opts);
  int nstates = state.capacity ();
  Matrix output (nsteps, nstates + 1);
  cout << "Computing Fitzhugh-Nagumo for times " << out_times(0) 
       << " to " << (int)out_times(out_times.length()-1) << endl; 
  output = ode.integrate (out_times);
  retval.resize (1);
  retval(0) = output;
  return retval;
}