数値解析の練習で，$\displaystyle f(x) = \frac{x^3}{e^x-1}$ や $g(x) = \dfrac{x^5}{e^x-1}$ の極大値を求める意義について。
単に演習のための演習問題というのではなく，以下のような意義があるんだよという話。
参考：
温度 $T$ の黒体から放射される放射強度は，周波数 $\nu$ の関数として単位周波数あたりであらわすと
\begin{eqnarray}
I(\nu) &=& \frac{2 h \nu^3}{c^2} \frac{1}{\exp\left(\frac{h\nu}{k T} \right) -1} \\
&=& \frac{2h}{c^2} \left( \frac{kT}{h}\right)^3 \left( \frac{h \nu}{kT} \right)^3 \frac{1}{\exp\left(\frac{h\nu}{k T} \right) -1} \\
&=& \mbox{const.} \times \frac{x^3}{e^x-1}, \quad x \equiv \frac{h\nu}{kT}
\end{eqnarray}
と書ける。ここで $h$ はプランク定数，$k$ はボルツマン定数，$c$ は光速。これがプランクの法則であり，プランクの法則であらわされる放射強度分布のことをプランク分布という。
つまり，関数 $f(x) \equiv \dfrac{x^3}{e^x-1}$ のグラフを描き，極大値を調べるということはプランク分布を知るという意義があるんだよ。
また，プランク分布を単位波長あたりにして波長 $\lambda$ の関数としてあらわすと
\begin{eqnarray}
I^{\prime}(\lambda) &=& \frac{2hc^2}{\lambda^5} \frac{1}{\exp\left(\frac{hc}{kT\lambda} \right)-1}\\
&=& 2 h c^2 \left( \frac{kT}{hc}\right)^5 \left( \frac{hc}{kT \lambda}\right)^5\frac{1}{\exp\left(\frac{hc}{kT\lambda} \right)-1}\\
&=& \mbox{const.} \times \frac{x^5}{e^x-1}, \quad x \equiv \frac{hc}{kT \lambda}
\end{eqnarray}
と書ける。
プランク分布 $I^{\prime}(\lambda)$ が最大となる波長 $\lambda_{\rm max}$ は温度 $T$ に反比例することが知られていて，これをウィーンの変位則という。
$I^{\prime}(\lambda)$ の最大値を調べることは $g(x) \equiv \dfrac{x^5}{e^x-1}$ の極大値を調べることだから，ウィーンの変位則を導くためにはこの関数が極大となる点を数値的に求める必要があるわけだ。
$g(x)$ が極大となる $x$ を $x_{\rm max}$ とすると，$g^{\prime}(x_{\rm max}) = 0$ を数値的に解き，
\begin{eqnarray}
x_{\rm max} &=& \frac{hc}{kT \lambda_{\rm max}} \\
\therefore\ \ \lambda_{\rm max} &=& \frac{hc}{k\, x_{\rm max}} \frac{1}{T} \propto \frac{1}{T}
\end{eqnarray}
となる。比例定数の具体的な数値も簡単に求められるよね。
ウィーンの変位則
$g(x) \equiv \dfrac{x^5}{e^x-1}$ が極大値をとる $x$ の値を $x_{\rm max}$ とすると，プランク分布が最大値をとる波長 $\lambda_{\rm max}$ は
\begin{eqnarray}
x_{\rm max} &=& \frac{hc}{kT \lambda_{\rm max}} \\
\therefore\ \ \lambda_{\rm max} &=& \frac{hc}{k\, x_{\rm max}} \frac{1}{T} \propto \frac{1}{T}
\end{eqnarray}
つまり，温度 $T$ に反比例する。これをウィーンの変位則という。ここで $h$ はプランク定数，$k$ はボルツマン定数，$c$ は光速。これら3つの定数はいずれも定義定数であり，誤差はない。
$x_{\rm max}$ を求める
ライブラリの import
# SymPy を使うときのおまじない
from sympy.abc import *
from sympy import *
# グラフは SymPy Plotting Backends (SPB) で描く
from spb import *
# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
import matplotlib.pyplot as plt
# mathtext font の設定
plt.rcParams['mathtext.fontset'] = 'cm'
$g(x)$ の定義
def g(x):
return x**5/(exp(x) - 1)
$g'(x)$
# g(x) の微分
dg = diff(g(x), x)
display(dg)
$g'(x)$ のグラフを描き，$g'(x) = 0$ となる $x$ のあたりをつける
xaxis = hline(0, rendering_kw={'c':'k'})
graphics(
line(dg, (x, 0.001, 10), r'$g^{\prime}(x)$'),
xaxis
);
上のグラフからわかるように，$g^{\prime}(x_{\rm max}) = 0$ となる $x_{\rm max}$ は $4 < x_{\rm max} < 6$ の範囲にある。
$g'(x_{\rm max}) = 0$ となる $x_{\rm max}$ を数値的に求める
dg = 0 となる $x$ を $4 < x < 6$ の範囲で数値的に求める。
x_max = nsolve(dg, x, [4, 6])
ウィーンの変位則の比例定数を求める
$$ W_{\rm const} \equiv \frac{h c}{k\, x_{\rm max}}$$
h = 6.62607015E-34
k = 1.380649E-23
c = 299792458
W_const = h * c/(k * x_max)
print('%e' % W_const)
太陽放射の $\lambda_{\rm max}$
太陽からの放射を表面温度 $T = 5800 \,\mbox{K}$ のプランク分布とすると，
$$ \lambda_{\rm max} = \frac{W_{\rm const}}{T}$$
lambda_max = W_const/5800
# 1 nm = 10^-9 m なので，ナノメートルにする。
nm = 1e-9
print('%d nm' % (lambda_max/nm))
この波長が，まさに人間の目で見える可視光の範囲に入っている事実に刮目せよ。