python中ODE的求解系统;在此调用上完成了多余的工作

问题描述

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

  div(class='sm:grid sm:grid-cols-3 sm:gap-4 sm:items-start sm:border-t sm:border-gray-200 sm:pt-5')
    label.block.text-sm.font-medium.text-gray-700(for='userPhotos' class='sm:mt-px sm:pt-2')
      | Review Photos
    .mt-1(class='sm:mt-0 sm:col-span-2')
      .max-w-lg.flex.justify-center.px-6.pt-5.pb-6.border-2.border-gray-300.border-dashed.rounded-md
        .space-y-1.text-center
          svg.mx-auto.h-12.w-12.text-gray-400(stroke='currentColor' fill='none' viewBox='0 0 48 48' aria-hidden='true')
            path(d='M28 8H12a4 4 0 00-4 4v20m32-12v8m0 0v8a4 4 0 01-4 4H12a4 4 0 01-4-4v-4m32-4l-3.172-3.172a4 4 0 00-5.656 0L28 28M8 32l9.172-9.172a4 4 0 015.656 0L28 28m0 0l4 4m4-24h8m-4-4v8m-12 4h.02' stroke-width='2' stroke-linecap='round' stroke-linejoin='round')
          .flex.text-sm.text-gray-600
            label.relative.cursor-pointer.bg-white.rounded-md.font-medium.text-indigo-600(for='userPhotos' class='hover:text-indigo-500 focus-within:outline-none focus-within:ring-2 focus-within:ring-offset-2 focus-within:ring-indigo-500')
              span Upload files
              input#userPhotos.sr-only(type='file' multiple,accept='image/*',name='userPhotos')
            p.pl-1 or drag and drop
          p.text-xs.text-gray-500
            | PNG,JPG,GIF up to 10MB

它不适用于(并打印 RuntimeWarning 消息):

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

知道第二段代码有什么问题吗?我对 python 比较陌生,不知道它的局限性。

解决方法

精华暗能量模型(见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]