2017-06-13 24 views
1

我有三個耦合的ODE,我在python中使用RK4方法來解決它們與給定的初始條件。當我運行它們時,出現以下錯誤:RuntimeWarning:exp中遇到的溢出

RuntimeWarning: overflow encountered in exp

對於我可能犯錯的地方的任何幫助將不勝感激。

下面的代碼(原貼here):

#simple exponential potential 
# u = K*phi'/H0; v = omega_matter**(1/3)*(1+z); w = l*K*phi' - ln((K**2)*V0/H0**2) 
# f,g,h are functions of derivation of u,v,w respectively derieved w.r.t t*H0 = T 

import matplotlib.pyplot as plt 
import numpy as np 
import math 

def f(T,u,v,w): 
    return -3*u*((v**3 + (u**2)/6 + np.exp(-w)/3)**0.5) + l*np.exp(-w) 

def g(T,u,v,w): 
    return -v*(v**3 + (u**2)/6 + np.exp(-w)/3)**0.5 

def h(T,u): 
    return l*u 

p = 0.1 
q = 1.0 
dh = 0.1 
n = (q-p)/dh 
u = [0.0] 
v = [1100] 
T = [0.00001] 
w = [-44.927] 
l = 1.3 

for i in range(0,int(n)): 
    k1 = f(T[i],u[i],v[i],w[i]) 
    r1 = g(T[i],u[i],v[i],w[i]) 
    print k1, r1 
    s1 = h(T[i],u[i]) 
    print s1 
    k2 = f(T[i] + 0.5*dh,u[i] + k1*dh,v[i] + k1*dh,w[i] + k1*dh) 
    r2 = g(T[i] + 0.5*dh,u[i] + r1*dh,v[i] + r1*dh,w[i] + r1*dh) 
    s2 = h(T[i] + 0.5*dh,u[i] + s1*dh) 
    print k2,r2,s2 
    k3 = f(T[i] + 0.5*dh,u[i] + k2*dh,v[i] + k2*dh,w[i] + k2*dh) 
    r3 = g(T[i] + 0.5*dh,u[i] + r2*dh,v[i] + r2*dh,w[i] + r2*dh) 
    s3 = h(T[i] + 0.5*dh,u[i] + s2*dh) 
    k4 = f(T[i] + dh,u[i] + dh*k3,v[i] + dh*k3,w[i] + k3*dh) 
    r4 = g(T[i] + dh,u[i] + r3*dh,v[i] + dh*r3,w[i] + r3*dh) 
    s4 = h(T[i] + dh,u[i] + dh*s3) 
    T == T.append(T[i] + dh) 
    u == u.append(u[i] + (dh/6)*(k1 + 2.0*k2 + 2.0*k3 + k4)) 
    v == v.append(v[i] + (dh/6)*(r1 + 2.0*r2 + 2.0*r3 + r4)) 
    w == w.append(w[i] + (dh/6)*(s1 + 2.0*s2 + 2.0*s3 + s4)) 

plt.plot(T,u, '-b') 
plt.plot(T,v, '-r') 
plt.plot(T,w, '-g') 
plt.title('quintessence cosmological model') 
plt.show() 
+1

歡迎您!嘗試用稍微不同的方式提出這個問題可能會有用。最好提供*最小*可能的代碼來重現相關問題,並在問題的主體中提供它。只是要求幫助調試大部分代碼的問題並不適合SO。您可以查看https://stackoverflow.com/help/how-to-ask,獲取有關如何在此提問的最佳指南。 – wphicks

+0

[overflow in exp,python]可能的重複(https://stackoverflow.com/questions/43100902/overflow-in-exp-python) –

+0

我強烈懷疑你有'f','g'和''錯了。要麼是這個問題,要麼是這個問題極其嚴重。在你的起點上'v'的導數值'-3619146406000.8389',已經看起來很大,經過一小步之後,你現在將巨大(負數)值送入'np.exp(-w)',並且結果導致無窮。 (我沒有看到RK4實現有什麼問題。)你有沒有參考你正在解決的物理問題? –

回答

1

你的問題似乎將在您的fg函數調用發起。當它被傳遞給該代碼中的numpy.exp函數時,請注意您的值爲-w。如果它的計算結果爲一個足夠大的正數,那麼調用exponential function似乎可能會導致溢出。谷歌告訴我,e^44.927將評估爲3.2474927e + 19。

您可能需要爲w選擇較小的絕對值才能運行此代碼而不會發生溢出錯誤。

(見Python-How to determine largest/smallest int/long/float/complex numbers my system can handle你如何能夠確定什麼最大允許的整數值是在Python)。

+0

int的最小和最大值在這裏並不真正相關;這都是浮點。 –

+0

@MarkDickinson好的。我提到了一個更好的問題。 –

+0

謝謝你的幫助。使用的代碼是正確的,但問題在於使用的微分方程。它們是自主的方程式,如果初始值很大,它會在一段時間後爆炸。關於如何計算自治方程式的任何想法或參考都會導致我無法在網絡中找到適當的解釋。 @MarkDickinson –