【发布时间】：2017-11-15 16:06:52
【问题描述】：
我有三个耦合的 ODE，我在 python 中使用 RK4 方法在给定的初始条件下求解它们。当我运行它们时，我收到以下错误：
RuntimeWarning: 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()
【问题讨论】：
-
欢迎您！尝试以稍微不同的方式提出这个问题可能会很有用。最好提供重现问题的最少 可能的代码，并在问题的正文中提供。仅在调试大部分代码时寻求帮助的问题并不适合 SO。您可以查看stackoverflow.com/help/how-to-ask，了解如何最好地在这里提问。
-
我强烈怀疑您的
f、g和h等式是错误的。要么，要么问题是极度病态的。v在您的起点处的导数具有值-3619146406000.8389，这看起来已经很大了，经过一小步之后，您现在将巨大的（负）值输入到np.exp(-w)中，结果是无穷大。 （我认为 RK4 的实现没有任何问题。）你有没有提到你正在解决的物理问题？ -
投票结束：这不是代码的问题，它可以很好地解决 ODE。这是方程式的问题。 （这可能是原始方程被翻译成 Python 的方式存在问题，但除非我们看到这些方程，否则我们无能为力。）
-
谢谢。下次会尝试跟随它。 @wphicks