[SciPy-user] Problems with scipy.integrate.ode

[Pearu Peterson]

On Sat, 19 Apr 2003, David Linke wrote:

> I am having problems with the set_f_params() method in ode.py.
> 
> An example showing the use of set_f_params() to supply additional 
> arguments (two or more) to f would be very heplful.
> 
> I tried a while to get it going but without success. I observed strange 
> behaviour when using a modified function test1 in ode.py to receive 
> additional arguments:
> 
> def test1():
>      #def f(t,y):
>      def f(t,y,b):
>          print type(b), b
>          a = Numeric.sin(6*t)
>          return y*y-a+y
>      ode_runner = ode(f)
>      ode_runner.set_integrator('vode')
>      ode_runner.set_initial_value([0.1,0.11,.1]*10)
                                                 ^^^ - ?!?
Note that ^^^ will return a list of 30 numbers.

>      par_b = 0.1,
>      ode_runner.set_f_params(par_b)
>      while ode_runner.successful() and ode_runner.t < 20:
>          y1 = ode_runner.integrate(ode_runner.t+2)
>          print ode_runner.t,y1[:3]
> 
> Is this the right way to use set_f_params?

Yes.

>  When I run this example the 
> arguments are passed to f only in the first two calls, as shown here:
> 
> Python 2.2.2 (#37, Oct 14 2002, 17:02:34) [MSC 32 bit (Intel)] on win32
>  >>> Integrators available: vode
> Found integrator vode
> <type 'tuple'> (0.10000000000000001,)        ! 1st call
> ...
> 2.0 [-0.1446766  -0.10226142 -0.1446766 ]
> <type 'tuple'> (0.10000000000000001,)        ! 2nd call
> ...
> 4.0 [-0.6351882 -0.5874518 -0.6351882]
> <type 'tuple'> (<NULL>,)                     ! 3rd call
> ...
> 6.0 [-0.93274299 -0.92029838 -0.93274299]
> <type 'tuple'> (<NULL>,)

I am also getting this in debian box. It looks like a ref-counting bug in
f2py callback hooks. I'll look for it.

Until then and also in general, you can pass additional parameters to
callback functions using optional arguments or global variables (this
works if there are no recursive calls to ode runner).

For example,

     par_b = [0.1]
     def f(t,y,b = par_b):
         print `b`
         a = Numeric.sin(6*t)
         return y*y-a+y
     ode_runner = ode(f)
     ode_runner.set_integrator('vode')
     ode_runner.set_initial_value([0.1,0.11,.1])
     while ode_runner.successful() and ode_runner.t < 10:
         y1 = ode_runner.integrate(ode_runner.t+2)
         print ode_runner.t,y1[:3]

I have used Python list for par_b because this will allow changing
`f` parameters during integration.

HTH,
	Pearu