[SciPy-User] multiple ODE problem with scipy.integrate.ode.ode

[Tony Yu]

On Mon, May 2, 2011 at 9:08 AM, Vitor Matos <vmatos at dei.uminho.pt> wrote:

> Hi,
> I've been trying to perform the integration of multiple ODEs.
>
> The problem is that only the first ODE being integrated produces the
> correct solution. The subsequent ODE solution come out as gibberish.
>
> To demonstrate the problem, I've written a simple code that solves two
> RC circuit ODEs:
>
> ----------------------------------------------------------------------
>
> from scipy.integrate import *
>
> def f(t, y, arg1):
>    # dV/dt = -1/RC (V - Vi)
>    return -1/(arg1[0]*arg1[1])*(y - arg1[2])
>
> y0 = 0
> t0 = 0
>
> # 0- R parameter
> # 1- C parameter
> # 2- Vin step
> arg = [100, 0.01, 5]
> arg2 = [200, 0.005, 5]
>
> r = ode(f).set_integrator('vode')
> r.set_initial_value(y0,t0)
> r.set_f_params(arg)
>
> r2 = ode(f).set_integrator('vode')
> r2.set_initial_value(y0,t0)
> r2.set_f_params(arg2)
>
> t_final = 10
> dt = 0.1
>
> t = [t0]
> v1 = [y0]
> v2 = [y0]
>
> while r.successful and r.t < t_final \
>        and r2.successful and r2.t < t_final:
>
>    r.integrate(r.t + dt)
>    r2.integrate(r2.t + dt)
>
>    #r2.integrate(r2.t + dt)
>    #r.integrate(r.t + dt)
>
>    print "1: ", r.t, r.y
>    print "2: ", r2.t, r2.y
>
>    t.append(r.t)
>    v1.append(r.y)
>    v2.append(r2.y)
>
> ## plot the simulation
> from matplotlib import *
> from pylab import *
> figure(1)    #Using matplotlib
> plot(t,v1,t,v2)
> title('RC circuit')
> xlabel('Time')
> show()
>
> ----------------------------------------------------------------------
>
> What could be the problem?
>
> Also, I'm new to Python and I'm still trying to grasp how to best use
> it.
>
> Thanks in advance,
> Vitor Matos
>

I'm not sure why this is happening, but if you separate the integration
loops (one for r followed by one for r2), the problem goes away. In other
words

while r.successful and r.t < t_final:
   r.integrate(r.t + dt)
   print "1: ", r.t, r.y
   t.append(r.t)
   v1.append(r.y)

while r2.successful and r2.t < t_final:
   r2.integrate(r2.t + dt)
   print "2: ", r2.t, r2.y
   t2.append(r2.t)
   v2.append(r2.y)

There's some weird sharing between the two calls to integrate (e.g. with a
single loop, moving r2.integrate before r.integrate results in r2
integrating smoothly, but r not). I don't know what's going on under the
hood of ode.integrate, so that's all the help I can give.

Best,
-Tony
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20110502/456d3d3b/attachment.html>