python 代码中耦合 ODE 的结果与 mathematica 不同

Question

根据我的语言知识，我的代码编写正确。但它没有给我正确的解决方案（情节）。当我在 mathematica 中解决相同的 ODE 系统时，我有正确的解决方案，两个解决方案完全不同。我正在写一个研究项目，所以我需要 python 中的适当代码。你能告诉我我的代码的错误吗？ python code solution Mathematica solution

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as si
##Three system
def func(state, T):
    H = state[0]
    P = state[1]
    R = state[2]
    Hd = -(16./3.)*np.pi*P 
    Pd = -4.*H*P
    Rd = H*R
    return Hd,Pd,Rd
T = np.linspace(0.1,0.9,50)
state0 = [1,0.0001, 0.1]
s = si.odeint(func, state0, T)
h = np.transpose(s)
plt.plot(T,h[0])
plt.show()

Mathematica 代码

Clear[H,\[Rho],a]
Eq1=(H'[t] == -16 \[Pi] \[Rho][t]/3)
Eq2= (\[Rho]'[t] == -4 H[t] \[Rho][t])
Eq3 = (a'[t] ==  H[t] a[t])
sol=NDSolve[{Eq1,Eq2, Eq3,
     H[0.1]==0.1, \[Rho][0.1]==0.1, a[0.1]==0.1}, 
       {H[t],\[Rho][t],a[t]}, {t,0.1, 0.9}]
Plot[Evaluate[{H[t]}/.sol],{t,0.1,0.9}]

Answer 1

两个代码都是正确的，我刚刚关闭我的笔记本电脑然后重新打开它，它给了我正确的结果（如 mathematica）

python 代码中耦合 ODE 的结果与 mathematica 不同

result of coupled ODE in python code is different from mathematica

python

wolfram-mathematica