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

[John W. Eaton]

On 23-Oct-2001, Douglas Eck <address@hidden> wrote:

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

Although it is nice to know how to do this (if you want to do it in a
stand-alone program, for example) I'm not sure that it is really
needed.  Most of what is slow about calling lsode from Octave is the
call to your RHS function when it is defined as an M-file.  So (I
think) you can get almost all of the speed increase if you just define
the RHS as a dynamically-linked function.  But maybe I'm missing
something.  Would anyone care to do some timings?

Thanks,

jwe


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



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------