计算科学 - 求解一组耦合 ODE 以获得正确的变量值 - 吾爱随笔录

求解一组耦合 ODE 以获得正确的变量值

计算科学 Python 颂 scipy 数值建模

2021-12-22 02:04:15

我的问题是关于如何解决 ODE 的耦合系统，并在图中打印出变量。

我正在求解一个 q 值和一个 e 值，在下面的这组耦合 ODE 中可以看到：

\begin{aligned} \frac{d q}{d t} & = \frac{48}{5 π M^{2}} (2 π M q)^{11 / 3} e \frac{1 + \frac{73}{24} e^{2} + \frac{37}{96} e^{4}}{(1 - e^{2})^{7 / 2}}, \\ \frac{d e}{d t} & = - \frac{304}{15 M} (2 π M q)^{8 / 3} e \frac{1 + \frac{121}{304} e^{2}}{(1 - e^{2})^{5 / 2}} . \end{aligned}

$\begin{aligned} \frac{dq}{dt} &= \frac{48}{5\pi M^2}(2\pi Mq)^{11/3}e\frac{1+\frac{73}{24}e^2+\frac{37}{96}e^4}{(1-e^2)^{7/2}},\\ \frac{de}{dt} &= -\frac{304}{15M}(2\pi Mq)^{8/3}e\frac{1+\frac{121}{304}e^2}{(1-e^2)^{5/2}}. \end{aligned}$

我的初始值为 $e$ 是 0.95 是我的初始值 $q$ 为 1，其中常数 $M$ = 1e9。

对于我的结果，我想绘制 q vs. e。我希望 e 值在 q 增加时变为零。

def q_e(x,t):
    # M, a constant
    M = 1e9

    # q and e input assignment
    q = x[0]
    e = x[1]

    # Eq.
    dqdt = (48/(5*math.pi*M**2))*(2*math.pi*M*q)**(11/3)*((1+(73/24)*e**2+(37/96)*e**4)/(1-e**2)**(7/2))
    dedt = -1*(304/(15*M))*(2*math.pi*M*q)**(8/3)*e*((1+(121/304)*e**2)/(1-e**2)**(5/2))
    return [dfdt,dedt]

x0 = [1,0.95]   # initial values q = 1, e = 0.99
t = np.linspace(0,100)
y = odeint(q_e,x0,t)

q = x[:,0]      # output in sep. columns 
e = x[:,1]

plt.plot(q,e)
#plt.semilogy(q,e) # test

我假设如果我打印出我的 q 和 e 值，我可以看到哪里出错了。于是我试了一下，

print(y)

并得到了输出

[[  1.00000000e+00   9.50000000e-01]
[  6.29120765e+05   3.44703571e-05]
[  1.39158341e+02   4.86785907e+02]
[  1.39158341e+02   4.86793360e+02]

其中 [ 1.39158341e+02 4.86793360e+02] 在之后的许多列中重复其长度为 $t$ . 注意到 [ 6.29120765e+05 3.44703571e-05] 的 q 和 e 值是正常的，这有点积极。

我试过弄乱我的代码，用 e 值代替时间 t，但没有任何效果。求指点和帮助，谢谢！

1个回答

功能 $q(e)$ 满足一阶线性 ODE

\frac{d q}{d e} = \frac{111 e^{4} + 876 e^{2} + 288}{(e^{2} - 1) (121 e^{2} + 304)} q (e),

$\frac{\mathrm{d}q}{\mathrm{d}e} = \frac{111 e^4+876 e^2+288}{(e^2-1) (121 e^2+304)} q(e),$ 这可以通过使用积分因子很容易地解决：

q (e) = C_{1} \exp (\frac{111}{121} e - 3 arctanh e + \frac{10440}{1331 \sqrt{19}} \arctan (\frac{11}{4 \sqrt{19}} e)) .

$q(e) = C_1 \exp\left( \tfrac{111}{121} e-3\operatorname{arctanh}e+\tfrac{10440}{1331 \sqrt{19}} \operatorname{arctan}\left(\tfrac{11}{4 \sqrt{19}} e\right) \right).$ 对于您给出的初始条件，这给出了

q |_{t = \infty} = q (0) / q (e_{0}) = 38.561177784533462867.

$q|_{t=\infty} = q(0)/q(e_0) = 38.561177784533462867.$

为了跟进我的评论，这就是重新缩放对数值解的作用（我添加了一个因子 $e$ 来dqdt匹配公式）。

from numpy import *
from scipy.integrate import odeint, solve_ivp

M = 1e9
scale = 1e22

def rhs(t, u):
  q, e = u[0], u[1]
  dqdt = (48/(5*math.pi*M**2))*(2*math.pi*M*q)**(11/3)*((1+(73/24)*e**2+(37/96)*e**4)/(1-e**2)**(7/2)) * e
  dedt = -1*(304/(15*M))*(2*math.pi*M*q)**(8/3)*e*((1+(121/304)*e**2)/(1-e**2)**(5/2))
  return [dqdt/scale, dedt/scale]

ic = [1.0, 0.95]

sol = solve_ivp(rhs, (0.0, 10.0), ic, method="LSODA", atol=0, rtol=1e-8)
print("success: ", sol.success)
print("final t: ", sol.t[-1])
print("final y: ", sol.y[:,-1])
print("y error: ", abs(sol.y[0,-1] - 38.561177784533462867))

通过重新缩放，它成功而没有警告。既然你说你有兴趣策划 $q(e)$ , 只需要增加 $t$ 直到 $e(t)$ 变得非常小，而且 $t$ 大大小于 $100$ . 输出：

success:  True
final t:  10.0
final y:  [3.85611807e+01 7.80753973e-19]
y error:  2.9253639937110165e-06

其它你可能感兴趣的问题