Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.python > #98375 > unrolled thread

Scipy odeint (LSODA) gives inaccurate results; same code fine in MATLAB ode15s/ode23s

Back to article view | Back to comp.lang.python

Contents

  Scipy odeint (LSODA) gives inaccurate results; same code fine in MATLAB ode15s/ode23s Abhishek <abhishek.mallela@gmail.com> - 2015-11-06 14:01 -0800
    Re: Scipy odeint (LSODA) gives inaccurate results; same code fine in MATLAB ode15s/ode23s Christian Gollwitzer <auriocus@gmx.de> - 2015-11-07 07:54 +0100
      Re: Scipy odeint (LSODA) gives inaccurate results; same code fine in MATLAB ode15s/ode23s Michael Torrie <torriem@gmail.com> - 2015-11-07 08:47 -0700

#98375 — Scipy odeint (LSODA) gives inaccurate results; same code fine in MATLAB ode15s/ode23s

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 generalized 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 have also specified small rtol and atol values, as well as a large mxstep. 

My code is below:

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

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

Thanks in advance for any help.

[toc] | [next] | [standalone]

#98387

Am 06.11.15 um 23:01 schrieb Abhishek:
> 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 generalized 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 have also specified small rtol
> and atol values, as well as a large mxstep.
>
> My code is below:
>
> [...Python code...]

>
> Why am I getting a "stiff-looking" graph in Python

What does "stiff-looking" mean? I only know stiff differential 
equations, but this leads to noisy results. Your plots are smooth.

> (http://i.stack.imgur.com/UGWSH.png), when MATLAB gives me a proper
> one (http://i.stack.imgur.com/F2jzd.jpg)? I would like to understand
> where the problem lies and how to solve it.

It is very hard to analyze such a problem, unless you also post the 
Matlab code and plot both solutions into a single graph.

What I can see at first is that the initial conditions can't be the 
same. Your Python graph starts around 10^7, while the Matlab graph 
starts at 10^5. WHat happens if you integrate a simpler system - say 
does it make sense to set some of the coefficients to zero and if you 
still get the same difference, can you simplify the program to expose 
the problem? How do you know that the Matlab code produces the correct 
answer?

	Christian

[toc] | [prev] | [next] | [standalone]

#98400

On 11/06/2015 11:54 PM, Christian Gollwitzer wrote:
> It is very hard to analyze such a problem, unless you also post the 
> Matlab code and plot both solutions into a single graph.

Also he may have a quicker response posting to the scipy list, where
scientists and mathematicians regularly use and talk about scipy:

http://www.scipy.org/scipylib/mailing-lists.html

[toc] | [prev] | [standalone]

Back to top | Article view | comp.lang.python

csiph-web