在同一组轴上多次绘制相同的函数,但参数不同 [Python]

计算科学 Python 绘图 效率
2021-12-28 09:46:27

我目前正在尝试绘制一个函数来描述宇宙学中不同世界模型的线性扰动增长。我希望能够在同一组轴上拥有所有曲线,但在设置它时遇到了困难。

这是功能

D(z)=5ΩMH022z1+zH3(z)dz

在哪里H(z)=H0[Ωm(1+z)3+ΩΛ]

我的目标是相对于 z 绘制这个函数 D,但是有多个密度参数不同的图(Ω)。

我管理了两个解决方案,但都不能完美运行,第一个效率很低(为每组参数添加新函数):

z = np.arange(0.0,10,0.1)

MOm = 1.0
MOv = 0.0

COm = 0.3
COv = 0.7

H0 = 70

def Mf(z):
    A = (5/2)*MOm*(H0**2)
    H = H0 * np.sqrt( MOm*((1+z)**3) + MOv )
    return A * ((1+z)/(H**3))

def MF(z):
    res = np.zeros_like(z)
    for i,val in enumerate(z):
        y,err = integrate.quad(Mf,val,np.inf)
        res[i] = y
    return res

def MD(z):
    return (H0 * np.sqrt( MOm*((1+z)**3) + MOv )) * MF(z)

def Cf(z):
    A = (5/2)*COm*(H0**2)
    H = H0 * np.sqrt( COm*((1+z)**3) + COv )
    return A * ((1+z)/(H**3))

def CF(z):
    res = np.zeros_like(z)
    for i,val in enumerate(z):
        y,err = integrate.quad(Cf,val,np.inf)
        res[i] = y
    return res

def CD(z):
    return (H0 * np.sqrt( COm*((1+z)**3) + COv )) * CF(z) 


plt.plot(z,MD(z),label="Matter Dominated")
plt.plot(z,CD(z),label="Current Epoch")

因此,我尝试使用 for 循环使其更简单,但无法弄清楚如何为循环内的每个绘图添加标签:

Om = (1.0,0.3)
Ov = (0.0,0.7)

for param1,param2 in zip(Om,Ov):
    def f(z):
        A = (5/2)*param1*(H0**2)
        H = H0 * np.sqrt( param1*((1+z)**3) + param2 )
        return A * ((1+z)/(H**3))

    def F(z):
        res = np.zeros_like(z)
        for i,val in enumerate(z):
            y,err = integrate.quad(f,val,np.inf)
            res[i] = y
        return res

    def D(z):
        return (H0 * np.sqrt( param1*((1+z)**3) + param2 )) * F(z)

    plt.plot(z,D(z))

有人可以帮助解释这样做的有效方法吗?或者如何使用 for 循环为绘图动态添加标签。任何帮助将不胜感激。

1个回答

以下代码在我的电脑上不到一秒的时间内给出了绘图。

from scipy.integrate import quad

def H(z, H0, MOm, MOv):
    return H0 * np.sqrt( MOm*((1+z)**3) + MOv )

# fixed values
COm = 0.3
COv = 0.7
H0 = 70

z_range = np.linspace(0, 1, 100)
Om_range = (1.0,0.3)
Ov_range = (0.0,0.7)

for (MOm, MOv) in zip(Om_range, Ov_range):
    d_range = [(5/2)*MOm*(H0**2) * 
               quad(lambda x: (1 + x)/H(z, H0, MOm, MOv)**3, z, np.inf)[0] for z in z_range]
    plt.plot(z_range, d_range, 
             label = '$\\Omega_m = ' + str(MOm)+ ', \\Omega_\\Lambda = ' + str(MOv)+'$')
plt.legend(loc = 0)
plt.xlabel('$z$')
plt.ylabel('$D(z)$')

在此处输入图像描述