我如何在数值上（有效地）积分和绘制我的函数？

Question

我想在下面绘制这个等式：

问题 1： 你看...因为我的函数是 ν 的函数，所以我必须计算域中每个 ν 的积分。我的问题是：最好的方法是什么？

我考虑过使用 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()

我如何在数值上（有效地）积分和绘制我的函数？

How can I numerically (and efficiently) integrate and plot my function?

python

numpy

graph

scipy

calculation