普朗克定律，频率数字答案

【问题标题】：Plancks Law, Frequency figures普朗克定律，频率数字
【发布时间】：2019-04-15 19:58:04
【问题描述】：

我想绘制普朗克定律的频率版本。我首先尝试独立完成此操作：

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
%matplotlib inline

# Planck's Law
# Constants
h = 6.62607015*(10**-34) # J*s
c = 299792458 # m * s
k = 1.38064852*(10**-23) # J/K
T = 20 # K
frequency_range = np.linspace(10**-19,10**19,1000000)

def plancks_law(nu):
    a = (2*h*nu**3) / (c**2)
    e_term = np.exp(h*nu/(k*T))
    brightness = a /(e_term - 1)
    return brightness

plt.plot(frequency_range,plancks_law(frequency_range))
plt.gca().set_xlim([1*10**-16 ,1*10**16 ])
plt.gca().invert_xaxis()

这不起作用，我在某种程度上遇到了缩放问题。我的下一个想法是尝试使用这个问题中的这个人的代码：Plancks Formula for Blackbody spectrum

import matplotlib.pyplot as plt
import numpy as np

h = 6.626e-34
c = 3.0e+8
k = 1.38e-23

def planck_f(freq, T):
    a = 2.0*h*(freq**3)
    b = h*freq/(k*T)
    intensity =  a/( (c**2 * (np.exp(b) - 1.0) ))
    return intensity

# generate x-axis in increments from 1nm to 3 micrometer in 1 nm     increments
# starting at 1 nm to avoid wav = 0, which would result in division  by zero.
wavelengths = np.arange(1e-9, 3e-6, 1e-9) 
frequencies = np.arange(3e14, 3e17, 1e14, dtype=np.float64) 

intensity4000 = planck_f(frequencies, 4000.)
plt.gca().invert_xaxis()

这不起作用，因为我得到了除以零的错误。除了我没有看到除以零之外，分母不应该是零，因为指数项不应该等于一。我选择了频率作为示例代码中波长值的转换。

谁能帮助解决这个问题或解释我如何才能得到普朗克定律的频率而不是波长？

【问题讨论】：

在第二个示例中，您似乎没有绘制值。这可能是为什么它被零除失败的原因吗？
如果我在第二个示例中添加plt.plot(frequencies, intensity4000)、plt.gca().set_xlim([3e14, 1e15]) 和plt.show()，它对我有用，尽管有一些溢出警告。正如@Asmus 正确指出的那样，警告表明您应该仔细查看数字的范围。
我想部分问题是frequency_range = np.linspace(10**-19,10**19,1000000) 应该是frequency_range = np.logspace(10**-19,10**19,1000000)。否则，您将计算10**13 点。我一般来说你需要在对数空间中操作（即，通过在计算中携带对数来找到强度对数），因为你会用大（或小）的指数和立方数来溢出你的浮点数。

标签： python numpy matplotlib

【解决方案1】：

您无法安全地处理如此大的数字；即使对于b = h*freq/(k*T) 的相对“小”值，您的float64 也会溢出，例如np.exp(709.)=8.218407461554972e+307 可以，但是np.exp(710.)=inf。您必须相应地调整单位（指数）以避免这种情况！

请注意，您链接到的the other question 也是如此，如果您在planck() 的定义中插入print( np.exp(b)[:10] )，您可以检查前十个评估的b，您会看到在前几次出现溢出。在任何情况下，只需在另一个问题中使用answer posted，但将plt.plot(wavelengths, intensity) 中的x 轴转换为频率（我希望你知道如何从一个到另一个）:-)

【讨论】：

但是，这些都不会导致第二个示例中的值被零除错误。他们更有可能只是忘记绘制它们或对未显示的值进行进一步处理。
@blubberdiblub 我想这可能取决于他们的 python/numpy 版本，我得到的是 RuntimeWarning: overflow encountered in exp intensity = a/( (c**2 * (np.exp(b) - 1.0) ))。
警告在默认设置下是正常的。但它会很高兴地继续使用无穷大作为值。默认设置下除以无穷大也可以，根据 a 的符号，分别产生 0.0 或 -0.0。
@blubberdiblub 抱歉，我应该换一种说法。在任何一种情况下，覆盖范围都会导致问题；在所讨论的第一个示例中，我们评估brightness = a /(e_term - 1) ，由于np.exp(0.0) 中的下溢，它产生除以零。在另一种情况下，由于np.exp() 变为inf，我们有一个溢出。
是的，你是对的，当然。你是对的，最好争取高质量的方程并且没有警告。不过，我的观点与此不同。在问题中，他们说 “这不起作用，因为我得到了除以零错误” 指的是第二个示例。但是，尽管有溢出警告，但我看不出它如何使用所提供的信息。受溢出影响的结果子集可能没用，也可能没用，但它没有失败。我想您可以将获得的 0 结果称为“不安全”，如果可以接受 0，则将其称为“安全”。