python中的ODE求解系统；在本次通话中完成的多余工作

Posted 2023-02-18

【中文标题】python中的ODE求解系统；在本次通话中完成的多余工作

【英文标题】：Solving system of ODEs in python; excess work done on this call

【发布时间】：2021-06-18 06:30:42

【问题描述】：

我正在尝试在 python 中解决一个耦合 ODE 系统以获得不同的潜力。它适用于特定类型的势（指数），但是一旦势由幂律描述，python 生成的图就完全不连贯，它经常只是为所有参数分配零值。我的编码适用于：

kr1 = 8*np.pi
#rho_m = a**(-3)
#V = np.e**(-kr1*x_1)
#dVdx = -kr1*np.e**(-kr1*x_1)

def quintessence (x, t):
    a = x[0]
    x_1 = x[1]
    x_2 = x[2]
    dadt = (kr1*a/np.sqrt(3))*np.sqrt(1/2 * x_2**2 + np.e**(-kr1*x_1) + a**(-3))
    dx_1dt = x_2
    dx_2dt = -np.sqrt(3)*kr1*np.sqrt(1/2 * x_2**2 + np.e**(-kr1*x_1) + a**(-3))*x_2 + kr1*np.e**(-kr1*x_1)

    return[dadt, dx_1dt, dx_2dt]

x0 = [0.0001, 0, 0]
t = np.linspace(0, 80, 1000)

x = odeint(quintessence, x0, t)

a = x[:,0]
x_1 = x[:,1]
x_2 = x[:,2]

plt.plot(t, a)
plt.show()

它不起作用（并打印 RuntimeWarning 消息）：

kr1 = 8*np.pi
#rho_m = a**(-3)
#V = M**2*x_1**(-2) with M=1
#dVdx = -2M**2*x_1**(-3) 

def quintessence2 (x, t):
    a = x[0]
    x_1 = x[1]
    x_2 = x[2]
    V = x_1**(-2)
    dVdx_1 = -2*x_1**(-3)
    dadt = (kr1*a/np.sqrt(3))*np.sqrt((1/2) * x_2**2 + V + a**(-3))
    dx_1dt = x_2
    dx_2dt = -np.sqrt(3)*kr1*np.sqrt((1/2) * x_2**2 + V + a**(-3))*x_2 + dVdx_1    
    return [dadt, dx_1dt, dx_2dt]

x0 = [.0001, 0.01, 0.01]
t = np.linspace(1, 80, 1000)

x = odeint(quintessence2, x0, t)

a = x[:,0]
x_1 = x[:,1]
x_2 = x[:,2]

plt.plot(t, a)
plt.show()

知道第二段编码可能有什么问题吗？我对 python 比较陌生，我不知道它的局限性。

【问题讨论】：

什么是 RuntimeWarning 消息？

请使用np.exp(u) 而不是np.e**u。您能否尝试使用去奇异化的三次方，而不是 u**(-3) 使用 u/(eps+u**4) 代替 eps=1e-4 或 1e-8？然后整合是否进行？

会不会是符号错误，最后一个词应该是...-dVdx_1？

@LutzLehmann 当我尝试去奇异化三次幂'ODEintWarning：要求的精度过高（公差太小）。以 full_output = 1 运行以获取定量信息。信息。和另一条消息“RuntimeWarning：在 sqrt del sys.path[0] 中遇到无效值”

@LutzLehmann 在原始等式中确实存在`-dVdx_1`，但正如您所见，当我定义 dVdx_1 时，我将减号放在那里。

【参考方案1】：

典型的暗能量模型（参见physics.SE作为任意参考）是

x'' + 3*H*x' + V'(x) = 0

H = a'/a

这里，根据所写的，

H = kr1/sqrt(3) * sqrt(E) with

E = 0.5*x'^2 + V(x) + a^(-3)

(( which leads to
   E' = x'*(x''+V'(x)) - 3*a^(-4)*a' = - 3*H*x'^2 -3*a^(-3)*H 
      = -3*H*(x'^2+a^(-3))
))

我会在这些部分实现这个系统，因为

def V(x): return ...
def dV(x): return ...

def quintessence(u,t):
    a,x,dx = u
    E = 0.5*dx**2 + V(x) + a**-3
    H = kr1/3**0.5 * E**0.5
    da = a*H
    ddx = -3*H*dx -dV(x)
    return [da, dx, ddx]