黑体辐射与蟒蛇的整合答案

【问题标题】：integration of Black Body radiation with python黑体辐射与蟒蛇的整合
【发布时间】：2021-02-08 04:42:14
【问题描述】：

我正在尝试绘制黑体辐射的积分。但由于某种原因，图表上有一个小飞艇。我也收到了警告。

IntegrationWarning：算法不收敛。舍入误差被检测到

我认为这可能高估了集成。但我不确定如何解决这个问题。

import scipy.integrate as integrate
import numpy as np                 #math module 

import matplotlib                  #Plotting module 
import matplotlib.pyplot as plt 
#constants
h = 6.626e-34
c = 3.0e+8
k = 1.38e-23           # Boltzmann constant
T = range(1, 6000)
plot = []

for t in T:
    integrand = lambda x : (((2*h*(c**2))/(x**5))*(1/(np.exp((h*c)/(x*k*t))-1)))
    result,err = integrate.quad(integrand, 0, np.inf, epsabs=0, limit=50)
    plot.append(result)
                
plt.figure(figsize=(10,10))        #Create a figure of a certain size
plt.plot(T,plot)                      #Make a plot in the figure.
plt.xlabel('Temperature K', fontsize=14)       #X Label 
plt.ylabel('Intensity (W/sr m^3)', fontsize=14)       #Y label
plt.title("Plank's Formula over all spectra")

【问题讨论】：

标签： python physics astronomy

【解决方案1】：

这是我的解决方法：

基本上，我不是从 0 积分到 np.inf，而是从小值积分到大值。我选择小值和大值 lambda_min 和 lambda_max 比参考 lambda 小/大 100 倍（使用 width），即 h * c / labmda = k_B * T。

另外，我绘制了黑体曲线的精确积分：sigma * T^4 作为一条线，它们彼此重叠。您可以增加width 以获得更准确的积分

import scipy.integrate as integrate
import numpy as np                 #math module 

import matplotlib                  #Plotting module 
import matplotlib.pyplot as plt 
plt.ion()
#constants
h = 6.626e-34
c = 3.0e+8
k = 1.38e-23           # Boltzmann constant
sigma = 5.6704e-8      #For later
width = 100            #For integrating. Increase for more accuracy
T = range(1, 6000)
plot = []

for t in T:
    lambda_rel = h * c / (k * t)     #defining the reference lambda
    lambda_min = lambda_rel / width  #Defining lower int bound
    lambda_max = lambda_rel * width  #Defining upper int bound
    integrand = lambda x : (((2*h*(c**2))/(x**5))*(1/(np.exp((h*c)/(x*k*t))-1)))
    result,err = integrate.quad(integrand, lambda_min, lambda_max, epsabs=0, limit=50)
    plot.append(result)
                
plt.figure(figsize=(10,10))        #Create a figure of a certain size
plt.plot(T,plot, 'r-')                      #Make a plot in the figure.
plt.plot(T, sigma * T**4, 'b-')    # Here I also plot the exact integral of the blackbody curve
plt.xlabel('Temperature K', fontsize=14)       #X Label 
plt.ylabel('Intensity (W/sr m^3)', fontsize=14)       #Y label
plt.title("Plank's Formula over all spectra")

【讨论】：