X-Received: by 10.66.102.65 with SMTP id fm1mr8760951pab.20.1446847269045; Fri, 06 Nov 2015 14:01:09 -0800 (PST)
X-Received: by 10.50.6.40 with SMTP id x8mr314103igx.0.1446847269004; Fri, 06 Nov 2015 14:01:09 -0800 (PST)
Path: csiph.com!au2pb.net!usenet.blueworldhosting.com!feeder01.blueworldhosting.com!peer01.iad.highwinds-media.com!news.highwinds-media.com!feed-me.highwinds-media.com!border1.nntp.dca1.giganews.com!nntp.giganews.com!i2no508411igv.0!news-out.google.com!fs1ni1058igb.0!nntp.google.com!i2no508394igv.0!postnews.google.com!glegroupsg2000goo.googlegroups.com!not-for-mail
Newsgroups: comp.lang.python
Date: Fri, 6 Nov 2015 14:01:08 -0800 (PST)
Complaints-To: groups-abuse@google.com
Injection-Info: glegroupsg2000goo.googlegroups.com; posting-host=129.237.108.11; posting-account=lO2z8woAAADplnjgj1WnWhHZZkr5y45P
NNTP-Posting-Host: 129.237.108.11
User-Agent: G2/1.0
MIME-Version: 1.0
Message-ID: <942d2bd0-6322-4733-a6dd-ceeed40368cb@googlegroups.com>
Subject: Scipy odeint (LSODA) gives inaccurate results; same code fine in MATLAB ode15s/ode23s
From: Abhishek <abhishek.mallela@gmail.com>
Injection-Date: Fri, 06 Nov 2015 22:01:09 +0000
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: quoted-printable
Lines: 137
X-Received-Bytes: 6037
X-Received-Body-CRC: 1361481398
Xref: csiph.com comp.lang.python:98375

I have code that runs perfectly well in MATLAB (using ode15s or ode23s) but=
 falters with Scipy odeint. The MATLAB code is for a specific case of the g=
eneralized Python code. Here I have tried to reproduce the specific case in=
 Python. The logic in the code is airtight and the algorithm is sound. I ha=
ve also specified small rtol and atol values, as well as a large mxstep.=20

My code is below:

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.integrate import odeint
   =20
    #Constants and parameters
    N =3D 2
    K00 =3D np.logspace(0,3,101,10)
    len1 =3D len(K00)
    epsilon =3D 0.01
    y0 =3D [0]*(3*N/2+3)
    u1 =3D 0
    u2 =3D 0
    u3 =3D 0
    Kplot =3D np.zeros((len1,1))
    Pplot =3D np.zeros((len1,1))
    S =3D [np.zeros((len1,1)) for kkkk in range(N/2+1)]
    KS =3D [np.zeros((len1,1)) for kkkk in range(N/2+1)]
    PS =3D [np.zeros((len1,1)) for kkkk in range(N/2+1)]
    Splot =3D [np.zeros((len1,1)) for kkkk in range(N/2+1)]
    KSplot =3D [np.zeros((len1,1)) for kkkk in range(N/2+1)]
    PSplot =3D [np.zeros((len1,1)) for kkkk in range(N/2+1)]
   =20
    for series in range(0,len1):
        K0 =3D K00[series]
        Q =3D 10
        r1 =3D 0.0001
        r2 =3D 0.001
        a =3D 0.001
        d =3D 0.001
        k =3D 0.999
        S10 =3D 1e5
        P0 =3D 1
        tfvec =3D np.tile(1e10,(1,5))
        tf =3D tfvec[0,0]
        time =3D np.linspace(0,tf,len1)
   =20
        #Defining dy/dt's
        def f(y,t):
            for alpha in range(0,(N/2+1)):
                S[alpha] =3D y[alpha]
            for beta in range((N/2)+1,N+1):
                KS[beta-N/2-1] =3D y[beta]
            for gamma in range(N+1,3*N/2+1):
                PS[gamma-N] =3D y[gamma]
            K =3D y[3*N/2+1]
            P =3D y[3*N/2+2]
   =20
            # The model equations
            ydot =3D np.zeros((3*N/2+3,1))
            B =3D range((N/2)+1,N+1)
            G =3D range(N+1,3*N/2+1)
            runsumPS =3D 0
            runsum1 =3D 0
            runsumKS =3D 0
            runsum2 =3D 0
   =20
            for m in range(0,N/2):
                runsumPS =3D runsumPS + PS[m+1]
                runsum1 =3D runsum1 + S[m+1]
                runsumKS =3D runsumKS + KS[m]
                runsum2 =3D runsum2 + S[m]
                ydot[B[m]] =3D a*K*S[m]-(d+k+r1)*KS[m]
   =20
            for i in range(0,N/2-1):
                ydot[G[i]] =3D a*P*S[i+1]-(d+k+r1)*PS[i+1]
   =20
            for p in range(1,N/2):
                ydot[p] =3D -S[p]*(r1+a*K+a*P)+k*KS[p-1]+ \
                          d*(PS[p]+KS[p])
   =20
            ydot[0] =3D Q-(r1+a*K)*S[0]+d*KS[0]+k*runsumPS
            ydot[N/2] =3D k*KS[N/2-1]-(r2+a*P)*S[N/2]+ \
                        d*PS[N/2]
            ydot[G[N/2-1]] =3D a*P*S[N/2]-(d+k+r2)*PS[N/2]
            ydot[3*N/2+1] =3D (d+k+r1)*runsumKS-a*K*runsum2
            ydot[3*N/2+2] =3D (d+k+r1)*(runsumPS-PS[N/2])- \
                            a*P*runsum1+(d+k+r2)*PS[N/2]
   =20
            ydot_new =3D []
            for j in range(0,3*N/2+3):
                ydot_new.extend(ydot[j])
            return ydot_new
           =20
        # Initial conditions
        y0[0] =3D S10
        for i in range(1,3*N/2+1):
            y0[i] =3D 0
        y0[3*N/2+1] =3D K0
        y0[3*N/2+2] =3D P0
   =20
        # Solve the DEs
        soln =3D odeint(f,y0,time,rtol =3D 1e-12, atol =3D 1e-14, mxstep =
=3D 50000000)
        for alpha in range(0,(N/2+1)):
            S[alpha] =3D soln[:,alpha]
        for beta in range((N/2)+1,N+1):
            KS[beta-N/2-1] =3D soln[:,beta]
        for gamma in range(N+1,3*N/2+1):
            PS[gamma-N] =3D soln[:,gamma]
   =20
        for alpha in range(0,(N/2+1)):
            Splot[alpha][series] =3D soln[len1-1,alpha]
        for beta in range((N/2)+1,N+1):
            KSplot[beta-N/2-1][series] =3D soln[len1-1,beta]
        for gamma in range(N+1,3*N/2+1):
            PSplot[gamma-N][series] =3D soln[len1-1,gamma]
   =20
        for alpha in range(0,(N/2+1)):
            u1 =3D u1 + Splot[alpha]
        for beta in range((N/2)+1,N+1):
            u2 =3D u2 + KSplot[beta-N/2-1]
        for gamma in range(N+1,3*N/2+1):
            u3 =3D u3 + PSplot[gamma-N]
   =20
        K =3D soln[:,3*N/2+1]
        P =3D soln[:,3*N/2+2]
        Kplot[series] =3D soln[len1-1,3*N/2+1]
        Pplot[series] =3D soln[len1-1,3*N/2+2]
        utot =3D u1+u2+u3
   =20
    #Plot
    plt.plot(np.log10(K00),utot[:,0])
    plt.show()

Why am I getting a "stiff-looking" graph in Python (http://i.stack.imgur.co=
m/UGWSH.png), when MATLAB gives me a proper one (http://i.stack.imgur.com/F=
2jzd.jpg)? I would like to understand where the problem lies and how to sol=
ve it.

Thanks in advance for any help.