2015-11-07 139 views
2

這個問題是對我之前的問題的一個後續處理:Assistance, tips and guidelines for converting Matlab code to Python成功將複雜的Matlab代碼轉換爲Python後,如何運行代碼?

我已經手動轉換了Matlab代碼。我正在使用MAC OS並從終端運行Python。但是,如何運行下面的代碼,獲得N的某個值,其中N是偶數?我應該得到一個圖(由劇情代碼指定)。

當我按原樣運行時,我什麼都沒有。

我的代碼如下:

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.integrate import odeint 

def Array(N): 
    K00 = np.logspace(0,3,101,10) 
    len1 = len(K00) 
    y0 = [0]*(3*N/2+3) 
    S = [np.zeros((len1,1)) for kkkk in range(N/2+1)] 
    KS = [np.zeros((len1,1)) for kkkk in range(N/2)] 
    PS = [np.zeros((len1,1)) for kkkk in range(N/2)] 
    Splot = [np.zeros((len1,1)) for kkkk in range(N/2+1)] 
    KSplot = [np.zeros((len1,1)) for kkkk in range(N/2)] 
    PSplot = [np.zeros((len1,1)) for kkkk in range(N/2)] 
    Kplot = np.zeros((len1,1)) 
    Pplot = np.zeros((len1,1)) 

    for series in range(0,len1): 
     K0 = K00[series] 
     Q = 10 
     r1 = 0.0001 
     r2 = 0.001 
     d = 0.001 
     a = 0.001 
     k = 0.999 
     P0 = 1 
     S10 = 1e5 
     tf = 1e10 
     time = np.linspace(0,tf,len1) 
     y0[0] = S10 
     y0[3*N/2+1] = K0 
     y0[3*N/2+2] = P0 
     for i in range(1,3*N/2+1): 
      y0[i] = 0 

     [t,y] = odeint(EqnsArray,y0,time, mxstep = 5000) 

     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-1] = y[:,gamma] 

     for alpha in range(0,(N/2+1)): 
      Splot[alpha][series] = y[len1-1,alpha] 
     for beta in range((N/2)+1,N+1): 
      KSplot[beta-N/2-1][series] = y[len1-1,beta] 
     for gamma in range(N+1,3*N/2+1): 
      PSplot[gamma-N-1][series] = y[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-1] 

     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) 
    plt.show() 

def EqnsArray(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-1] = 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] 
     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] 

    for p in range(1,N/2): 
     ydot[p] = -S[p]*(r1+a*K+a*P)+k*KS[p-1]+d*(PS[p-1]+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-1] 
    ydot[G[N/2-1]] = a*P*S[N/2]-(d+k+r2)*PS[N/2-1] 
    ydot[3*N/2+1] = (d+k+r1)*runsumKS-a*K*runsum2 
    ydot[3*N/2+2] = (d+k+r1)*(runsumPS-PS[N/2-1])- \ 
        a*P*runsum1+(d+k+r2)*PS[N/2-1] 

    ydot_new = [] 
    for j in range(0,3*N/2+3): 
     ydot_new.extend(ydot[j]) 
    return ydot_new 
+2

你必須在你的代碼最後調用你的函數。 – jmugz3

回答

3

你打電話給你的功能,如:

Array(12) 

你必須在你的代碼的末尾添加此。

+0

具體而言,將此行添加到程序的最後。 – Adam

+0

好吧,我得到一個「全球名稱」N'未定義「錯誤,指def EqnsArray(y,t)後的行。我會在哪裏放置全局命令? – abscissa

+0

您需要將'N'作爲參數添加到您的'EqnsArray'函數中。 – dimo414