[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Calling lsode directly from a c++ program via liboctave?
From: 
John W. Eaton 
Subject: 
Re: Calling lsode directly from a c++ program via liboctave? 
Date: 
Thu, 8 Nov 2001 23:37:03 0600 
On 23Oct2001, 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 dynamicallylinked 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
standalone 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 Mfile. So (I
think) you can get almost all of the speed increase if you just define
the RHS as a dynamicallylinked 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
 FitzhughNagumo 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 twovarible
 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
 FitzhughNagumo 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 twovarible
 simplification of the famous Hodgkin Huxley equations.

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

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

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

 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/ovstruct.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,,
 "FitzhughNagumo 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 FitzhughNagumo 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
