数值求解 Friedmann-Lemâitre 方程时的问题答案

【问题标题】：Problem when numerically solving Friedmann-Lemâitre equation数值求解 Friedmann-Lemâitre 方程时的问题
【发布时间】：2021-08-30 19:09:41
【问题描述】：

我正在尝试求解 a=1 到 a = 0 的 FL 方程。但是，当我尝试使用 odeint 求解时，在 a 和 t 的值较小时，我的平方根中会出现负值衍生品，把事情搞砸了。是我使用了错误的方法来解决它还是我的方程式中遗漏了一些东西？

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

#FL equation
#(da/dt/a)^2 = H0^2*[(8.4*10^-5)/a^4 + 0.31/a^3  +  0.69 + (sum(omegas)-1)/a^2 )
H0=7.2e-11
omegaR = 8.4*10**-5
omegaM = 0.31
omegaA = 0.69

def model(a,t):
    RHS = H0**2*(omegaR/a**2 + omegaM/a + omegaA*a**2 +(omegaM+omegaA+omegaR-1) )
    dadt = np.sqrt(RHS)
    return dadt

# initial condition
y0 = 1

# time points
t = np.linspace(0,-10**14,2*10**7)

# solve ODE
a = odeint(model,y0,t)

At low values of t, it just cuts off, which should not happen

【问题讨论】：

请提供预期的minimal, reproducible example (MRE)。我们应该能够复制和粘贴您的代码的连续块，执行该文件，并重现您的问题以及跟踪问题点的输出。这让我们可以根据您的测试数据和所需的输出来测试我们的建议。显示中间结果与您的预期不同的地方。
您发布的代码无法运行：您缺少所有导入，并且未能定义odeint。
我只是忘了粘贴它们，我将问题编辑为包含所有内容并删除了非必要代码。

标签： python numerical-integration

【解决方案1】：

我相信你会被计算错误杀死。首先，您仍然没有采纳我的建议来追踪您的问题。 model 中的一个简单的print，就在故障点之前，显示了当时出了什么问题：

print(RHS, a**2, omegaR, omegaM, omegaA, omegaM+omegaA+omegaR-1)

这会得到一个输出流，显示 a 缓慢减少，直到我们最终得到：

[3.21061361e-21] [7.33441396e-08] 8.400000000000001e-05 0.31 0.69 8.399999999997299e-05
[2.08853856e-22] [7.34183571e-08] 8.400000000000001e-05 0.31 0.69 8.399999999997299e-05
[-2.25330684e-21] [7.34793736e-08] 8.400000000000001e-05 0.31 0.69 8.399999999997299e-05

您的 RHS 数字非常接近 0，以至于数值大小的巨大差异以及二进制表示中的微小错误都会给您带来负面结果。

更改您的操作顺序以首先保留最低有效数字。

【讨论】：