我写了以下代码。它是一个ODE,其中包含一个参数作为另一个ODE。正如我们所看到的,M(m0,z,b,c)在另一个ODE中使用,它本身就是一个ODE函数。代码很慢,有人能给我一个改进的建议吗?
import numpy as np
from scipy.integrate import odeint
def model(m,z,c,b):
dmdt = ((c**2-m)/(1+z))*(6-9*(m/c**2)+3*b*(m+(m**2)))
return dmdt
def M(m0,z,c,b):
m = odeint(model,m0,[0,z], args= (c, b))
mm=m[-1,0]
return mm
def model1(H ,z,m0,c,b):
c = 0.6
b=0.035
dHdt = (H/(1+z))*(6-9*(M(m0,z,c,b)/c**2)+3*b*(M(m0,z,c,b)+(M(m0,z,c,b)**2)))
return dHdt
def model2(H0,z,m0,c,b):
H = odeint(model1,H0,[0,z], args=(m0,c,b))
HH=H[-1,0]
return HH
print(model2(70,1,0.75,0.69,0.035))发布于 2018-05-02 16:26:15
您可以将耦合系统作为耦合系统来求解。
def model(U,z,c,b):
M, H = U
dMdt = ((c**2-M)/(1+z))*(6-9*(M/c**2)+3*b*(M+M**2))
dHdt = (H /(1+z))*(6-9*(M/c**2)+3*b*(M+M**2))
return [dMdt, dHdt]
def solution(H0,z,m0,c,b):
U = odeint(model,[m0,H0],[0,z], args=(c,b))[-1]
M, H = U
return H
print(solution(70,1,0.75,0.69,0.035)),它在代码修改后迅速返回0.107569653042。
def model1(H, z, m0, c, b):
mm = M(m0,z,c,b)
dHdt = (H/(1+z))*(6-9*(mm/c**2)+3*b*(mm+(mm)**2)))
return dHdt返回类似的0.107569746892,稍微慢一点。这6位符合的数字与1e-6的默认误差公差一致。
为了获得精度较高的结果,对误差容差atol, rtol的控制参数进行了设置。
为进一步减少业务,请
def model(U,z,c,b):
M, H = U
factor = (6-9*M/c**2+3*b*(M+M**2))/(1+z)
return [(c**2-M)*factor, H*factor]如果您的任务非常庞大,请使用编译后的编程语言进行快速批量数字处理。
https://stackoverflow.com/questions/50136298
复制相似问题