python中ODE的求解系统；在此调用上完成了多余的工作 [英] Solving system of ODEs in python; excess work done on this call

问题描述

我正在尝试在 python 中解决不同电位的耦合 ODE 系统.它适用于特定类型的势(指数)，但是一旦势由幂律描述，python 生成的图就完全不连贯，它经常只为所有参数分配零值.我的编码适用于:

I am trying to solve a system of coupled ODEs in python for different potentials. It works for a particular type of potential (exponential) but once the potential is described by a power law, the graph produced by python is not at all coherent and it frequently just assigns zero value to all arguments. My coding works for:

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 消息):

It doesnt work for (and prints a RuntimeWarning message):

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 比较陌生，不知道它的局限性.

Any idea what may be wrong with the second piece of coding? I am relatively new to python and I don't know its limitations.

推荐答案

quintessence 暗能量模型(参见 physics.SE 作为任意参考)是

The quintessence dark energy model (see physics.SE as an arbitrary reference) is

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))
))

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

I would implement this system in exactly these parts, as

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]

这篇关于python中ODE的求解系统；在此调用上完成了多余的工作的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！