2013-10-21 44 views
0

我有兩個程序,一個可以採用N個耦合ODE,另一個採用2個耦合ODE。如果我在兩個代碼中輸入2個相同的ODE,並且時間跨度相同,我會得到不同的答案。我知道正確的答案,所以我可以推斷出我的許多程序是錯誤的。爲什麼odeint時間在python中演變這兩種不同?

這裏是專用一個用於2方程的代碼:

# solve the coupled system dy/dt = f(y, t) 
def f(y, t): 
    """Returns the collections of first-order 
    coupled differential equations""" 
    #v11i = y[0] 
    #v22i = y[1] 
    #v12i = y[2] 
    print y[0] 

    # the model equations 
    f0 = dHs(tRel,vij)[0].subs(v12,y[2]) 
    f1 = dHs(tRel,vij)[3].subs(v12,y[2]) 
    f2 = dHs(tRel,vij)[1].expand().subs([(v11,y[0]),(v22,y[1]),(v12,y[2])]) 

    return [f0, f1, f2] 


# Initial conditions for graphing 
v110 = 6    
v220 = 6     
v120 = 4     
y0 = [v110, v220, v120]  # initial condition vector 

sMesh = np.linspace(0, 1, 10e3) # time grid 

# Solve the DE's 
soln = odeint(f, y0, sMesh) 

和這裏在N方程專用之一:

def f(y, t): 
    """Returns the derivative of H_s with initial 
    values plugged in""" 
    # the model equations 
    print y[0] 

    for i in range (0,len(dh)): 
     for j in range (0,len(y)): 
      dh[i] = dh[i].subs(v[j],y[j]) 

    dhArray = [] 
    for i in range(0,len(dh)): 
      dhArray.append(dh[i]) 

    return dhArray 

sMesh = np.linspace(0, 1, 10e3) # time grid 
dh = dHsFunction(t, V_s).expand() 
soln = odeint(f, v0, sMesh) 

其中DHS(TREL,維吉)= dHsFunction(T, V_s)即完全相同的ODE。同樣,y0和v0完全相同。但是,當我打印Y [0]在N許多情況下,我得到的輸出:

6.0 
5.99999765602 
5.99999531204 
5.97655553477 
5.95311575749 
5.92967598021 
5.69527820744 
5.46088043467 
5.2264826619 
2.88250493418 
0.53852720647 
-1.80545052124 
-25.2452277984 
-48.6850050755 
-72.1247823527 
-306.522555124 

而非的2專用的情況下:

6.0 
5.99999765602 
5.99999765602 
5.99999531205 
5.99999531205 
5.98848712729 
5.98848712125 
5.97702879748 
5.97702878476 
5.96562028875 
5.96562027486 
5.91961750442 
5.91961733611 
5.93039037809 
5.93039029335 
5.89564277275 
5.89564273736 
5.86137647436 
5.86137638807 
5.82758984835 

etc. 

其中所述第二結果是正確的一個和繪製正確的圖。

請讓我知道是否需要更多的代碼或其他。謝謝。

回答

1

您的第二個版本f會修改全局變量dh的值。 在第一次調用時,將替換其中的值,然後在隨後的所有調用中使用這些相同的值。

避免使用例如dh_tmp = list(dh)裏面的函數。

+0

我幾乎可以肯定你是對的。我試圖遍歷dh的每個元素,以便用變量替換它的值,如果它存在的話。建議採用其他方式來做到這一點。我試圖設置一個tmp = dh,並用tmp替換所有的dh,並得到相同的結果。 – faceforest

+0

閱讀http://docs.python.org/2/tutorial/classes.html#a-word-about-names-and-objects可能會澄清問題。 –

相關問題