CODEWITHSUNDEEP
pythonscipydifferential-equations
Francesco Di Lauro
Francesco Di Lauro

Reputation: 630

python odeint gives strange results

I am trying to understand how scipy.odeint works, but I have some problems. For exampe:

import numpy as np
from scipy.integrate import  odeint
from scipy.integrate import  odeint

def Diffeq(v,t, lam, gam,c, a):
    vdot = v
    for i in range(0,len(v)):
        if i == 0:
            vdot[0] = c[1]*v[1]- lam*a[0]*v[0]
        elif i == (len(v)-1):
            vdot[i] =  lam*a[i-1]*v[i-1] - (lam*a[i]+c[i])*v[i]
        else:
            vdot[i]=  lam*a[i-1]*v[i-1] - (lam*a[i]+c[i])*v[i]+ c[i+1]*v[i+1]
    print vdot
    return vdot



incond=np.array([0]*900)
incond[1] =1
t = np.linspace(0.0, 2, 1000)
ak = [2*i for i in range(0,900)]
lamma =2
gamma =1
c=[i*gamma for i in range(0,900)]

y = odeint(Diffeq, incond, t, args=(lamma,gamma,c,ak) )

This code should compute a system of differential equations of the form: enter image description here

where

xdot_0 = -(a_0 + c_0)*x_0(t) + c_1*x_1(t)
xdot_899 = a_898*x_898(t) -(a_899 + c_899)*x_900(t)

with initial condition x(0) = (0,1,0...0) when I try to analyse the results, I notice that my functions explode to + infinity. If I play with the constants ak, lamma, and gamma, I can get the results stuck in something like

x_k(t) = [0,-21,21, 0, 0,...,0]

at each operation. Therefore I am thinking that I made some mistake in my code, but I do not see where.

Upvotes: 0

Views: 361

Answers (1)

Warren Weckesser
Warren Weckesser

Reputation: 114781

In Python, when you execute the line

    vdot = v

vdot is not a copy of v. The two names now refer to the same object. So when you modify vdot in the function Diffeq, you also modify the input argument. If you modify vdot[0] and then later try to use v[0], you are actually getting the modified value vdot[0], so your calculations are not correct.

Change that line to, say,

    vdot = np.empty_like(v)

When I do that (and I remove the print statement, so the function finishes in a reasonable amount of time), odeint returns successfully. Here's a plot of the solution components:

plot

Upvotes: 2

Related Questions