python - 我如何在数值上(并有效地)整合和绘制我的函数？

Question

我想在下面绘制这个方程的图:

问题一:你看...因为我的函数是 ν 的函数，所以我必须计算我域中每个 ν 的积分。我的问题是:最好的方法是什么？
我想过使用 scipy 进行积分并使用 for 循环对每个 ν 计算几次，但这似乎是解决我的问题的一种非常不雅的方法。有人知道更好的选择吗？有人有不同的想法吗？
问题2:当我编写代码时，我遇到了一些错误，主要是因为我认为指数的指数非常小。您对我应该如何更改它以便我可以使用 Python 绘制此图有任何想法吗？
哦，如果你用不同的方法尝试，它应该看起来like this
这是我正在处理的代码。我现在又回到 Python 上了，所以可能有一些错误。我得到的情节与这应该看起来的情节非常不同。

from scipy.integrate import quad
from scipy.constants import c, Planck, k, pi
import numpy as np
import matplotlib.pyplot as plt

def luminosity_integral(r, x):
    T_est = 4000 
    R_est = 2.5 * (696.34*1e6)

    Temp = ((2/(3*pi))**(1/4)) * T_est * ((R_est/r)**(3/4))
    termo1 = ((4 * (pi**2) * Planck * (x**4) ) / (c**2))
    termo2 = ((Planck * x) / (k*Temp))
   

    return ((termo1 * r ) / (np.exp(termo2) - 1))


freqs = np.linspace(1e10, 1e16)
y = np.array([])
for i in freqs:
I = quad(luminosity_integral, (6 * 2.5 * (696.34*1e6)), (7e4 *  2.5 * (696.34*1e6)), args = (i))
temp = np.array([I[0]])
y = np.concatenate((y, temp))


plt.loglog(freqs, y)
plt.show()

Answer 1

重用术语 R_est 而不是将其表达式写入 3 次(如果您想更改该参数更好)。

您使用了 pi**2在常数乘以积分(不影响形状)

形状类似于您作为引用放置的形状，但不在建议的范围内。

您正在使用 T 的值如 T_* ， 你确定吗？

试试这个版本的代码

from scipy.integrate import quad
from scipy.constants import c, Planck, k, pi
import numpy as np
import matplotlib.pyplot as plt

R_est = 2.5 * (696.34e6)
def luminosity_integral(r, x):
    T_est = 4000 

    termo1 = ((4 * pi * Planck * (x**4) ) / (c**2))
    termo2 = ((Planck * x) / (k*T_est)) * (3*pi/2 * (r/R_est)**3)**0.25
    termo3 = np.exp(-termo2)
    return ((termo1 * r ) * termo3 / (1 - termo3))


freqs = np.logspace(6, 16)
y = np.zeros_like(freqs)

for i, nu in enumerate(freqs):
    y[i] = quad(luminosity_integral, (6* R_est), (7e4 *  R_est), args = (nu))[0]
    

plt.loglog(freqs, y)
plt.ylim([1e6, 1e25])
plt.show()