第43章波動方程式と熱伝導方程式 — 導出と変数分離法
第42章では,偏微分方程式の一般論と,ラグランジュの1階偏微分方程式の解き方を学んだ.本章からは,いよいよ自然界そのものを記述する具体的な偏微分方程式——波動方程式(wave equation)と熱伝導方程式(heat equation)——に取り組む.どちらも未知関数を2回・1回と偏微分した「2階線形偏微分方程式」であり,弦の振動,太鼓の膜,金属棒の温度分布,そして電磁波の伝播までを,たった1本の式で記述してしまう.本章の目標は2つある.(1) これらの方程式が物理法則(Maxwell方程式,運動方程式,フーリエの熱伝導法則)から実際にどう導かれるのかを,式変形を1行も飛ばさずに追うこと.(2) 導かれた方程式を,境界条件・初期条件のもとで実際に解く最強の道具——変数分離法(method of separation of variables)——を身につけることである.
変数分離法の考え方は単純明快である.「$x$ と $t$ の両方に依存する関数 $u(x,t)$ を,$x$ だけの関数 $X(x)$ と $t$ だけの関数 $T(t)$ の積 $u=X(x)T(t)$ という特別な形に限って探してみよう」というのである.なぜそんな特殊な形に限定してよいのかと思うかもしれないが,方程式が線形であるおかげで,見つかった特殊な形の解をいくらでも足し合わせた「無限1次結合」もまた解になる(重ね合わせの原理).そして,境界条件を満たす特殊な形の解(固有関数)の無限1次結合の係数を,フーリエ級数展開(第28章・第29章)によって決定すれば,任意の初期条件に対応できる一般的な解が完成する.この「変数分離 → 固有値問題 → フーリエ係数決定」という3段階の流れは,本章の4つの例題(1次元・2次元の波動方程式と熱伝導方程式)で繰り返し使われ,第45章以降のシュレーディンガー方程式でも全く同じ骨格が使われることになる.本章はいわば,量子力学の数学的土台を先取りする章でもある.
- 真空中のMaxwell方程式から,電場が満たす波動方程式 $\nabla^2\bm E=\dfrac{1}{c^2}\pdiff{^2\bm E}{t^2}$ を導く(ベクトル恒等式 $\nabla\times(\nabla\times\bm E)=-\nabla^2\bm E+\nabla(\nabla\cdot\bm E)$ を使う)
- (番外)波動方程式に $u=\psi(x)\cos(\omega t+\theta_0)$ を代入し,時間に依存しないシュレーディンガー方程式を導く——第45章への橋渡し
- 弦の微小要素の運動方程式 $ma=F$ から,1次元波動方程式 $\pdiff{^2u}{x^2}=\dfrac1{v^2}\pdiff{^2u}{t^2}$($v^2=T/\rho$)を導く
- 棒の微小区間の熱収支とフーリエの熱伝導法則から,熱伝導方程式 $\pdiff ut=a\pdiff{^2u}{x^2}$($a=k/(\rho C)$)を導く
- 変数分離法 $u=X(x)T(t)$ と,重ね合わせの原理,フーリエ正弦級数によるフーリエ係数の決定という3段階の解法パターン
- 1次元・2次元の波動方程式,1次元・2次元の熱伝導方程式を,具体的な境界値問題(両端固定,三角形分布・矩形分布の初期条件など)として実際に解き切る
- 波動方程式の解(振動し続ける)と熱伝導方程式の解(指数的に減衰する)の質的な違いを,同じ初期条件から出発して比較する
- 2次元ラプラス方程式 $\pdiff{^2u}{x^2}+\pdiff{^2u}{y^2}=0$ の境界値問題(ノートで中断している計算を本章で完成させる)
もとにしたノート:望月泰英『数学ノート 偏微分方程式』 pp. 11–22.
43.1 波動方程式の導出(Maxwell方程式から)
43.1.1 真空中のMaxwell方程式と電場の波動方程式
電気と磁気に関するすべての現象は,4本のMaxwell方程式(Maxwell's equations)にまとめられている.電場を $\bm E$,電束密度を $\bm D$(多くの物質では $\bm D=\varepsilon\bm E$ で,$\varepsilon$ は誘電率),磁束密度を $\bm B$,磁界の強さを $\bm H$(多くの物質では $\bm B=\mu\bm H$ で,$\mu$ は透磁率),電荷密度を $\rho$,電流密度を $\bm j$ と書くと,4本の方程式は
$$ \nabla\cdot\bm D=\rho, \qquad \nabla\cdot\bm B=0, \qquad \nabla\times\bm E=-\pdiff{\bm B}{t}, \qquad \nabla\times\bm H=\bm j+\pdiff{\bm D}{t} $$である.$\nabla\cdot(\cdot)$(発散,divergence)と $\nabla\times(\cdot)$(回転,rotation または curl)の定義は第19章で学んだ.物理的な意味を一言ずつ確認しておこう.第1式(Gaussの法則)は「電荷があるところから電気力線が湧き出す」こと,第2式は「磁気力線には湧き出し口(=孤立した磁極,モノポール)が存在しない」こと,第3式(Faradayの電磁誘導の法則)は「磁束密度が時間変化すると,それを取り囲むように電場が渦を巻いて生じる」こと,第4式(Ampère–Maxwellの法則)は「電流,または電束密度の時間変化(変位電流)が,それを取り囲むように磁場の渦を作る」ことを表している.
本節では,何もない真空中(電荷も電流も存在しない領域)を考える.つまり $\rho=0$,$\bm j=0$ とおく.真空中では $\bm D=\varepsilon_0\bm E$,$\bm B=\mu_0\bm H$($\varepsilon_0$:真空の誘電率,$\mu_0$:真空の透磁率)という比例関係(構成方程式)が成り立つので,4本の方程式は次のように書き換えられる.
$$ \nabla\cdot\bm E=0, \qquad \nabla\cdot\bm B=0, \qquad \nabla\times\bm E=-\pdiff{\bm B}{t}, \qquad \nabla\times\bm H=\pdiff{\bm D}{t} $$目標は,この4本の連立方程式から $\bm H$(や $\bm B$)を消去して,電場 $\bm E$ だけが満たす方程式を作ることである.鍵になるのは,任意のベクトル場 $\bm E$ について成り立つ次のベクトル恒等式である(第19章の発散・回転の性質から導かれる恒等式で,証明は第19章の演習または付録を参照).
\begin{equation} \nabla\times(\nabla\times\bm E)=-\nabla^2\bm E+\nabla(\nabla\cdot\bm E) \label{eq:43-vecid} \end{equation}ここで $\nabla^2$ はラプラシアン(Laplacian,$\nabla^2 f=\pdiff{^2f}{x^2}+\pdiff{^2f}{y^2}+\pdiff{^2f}{z^2}$,ベクトル場に対しては各成分に作用する)であり,$\nabla(\nabla\cdot\bm E)$ は「$\bm E$ の発散(スカラー量)」の勾配(gradient)である.
導出:電場の波動方程式
真空中では $\nabla\cdot\bm E=0$ であったから,恒等式\eqref{eq:43-vecid}の右辺第2項は消えて,
\begin{equation} \nabla\times(\nabla\times\bm E)=-\nabla^2\bm E \label{eq:43-step1} \end{equation}となる.一方,左辺の $\nabla\times\bm E$ には第3式(Faradayの法則)が使えて,$\nabla\times\bm E=-\pdiff{\bm B}t=-\mu_0\pdiff{\bm H}t$($\bm B=\mu_0\bm H$ を代入した)である.したがって左辺は
$$ \nabla\times(\nabla\times\bm E)=\nabla\times\left(-\mu_0\pdiff{\bm H}{t}\right) $$と書き直せる.ここで,$\mu_0$ は定数だから $\nabla\times$(空間についての微分)の外に出せて,さらに空間微分 $\nabla\times$ と時間微分 $\pdiff{}{t}$ は独立な変数についての微分なので順序を交換できる(これは高校で学ぶ「偏微分の順序交換」と同じ発想であり,第6章で学んだ性質である).よって
$$ \nabla\times\left(-\mu_0\pdiff{\bm H}{t}\right)=-\mu_0\pdiff{}{t}(\nabla\times\bm H) $$となる.ここに第4式 $\nabla\times\bm H=\pdiff{\bm D}t=\varepsilon_0\pdiff{\bm E}t$(真空中の構成方程式 $\bm D=\varepsilon_0\bm E$ を代入した)を代入すると,
$$ -\mu_0\pdiff{}{t}(\nabla\times\bm H)=-\mu_0\pdiff{}{t}\left(\varepsilon_0\pdiff{\bm E}t\right)=-\varepsilon_0\mu_0\pdiff{^2\bm E}{t^2} $$が得られる.以上をまとめると,式\eqref{eq:43-step1}の左辺は $-\varepsilon_0\mu_0\pdiff{^2\bm E}{t^2}$ に等しいことが分かったので,
$$ -\nabla^2\bm E=-\varepsilon_0\mu_0\pdiff{^2\bm E}{t^2} $$両辺に $-1$ を掛けて整理すると
$$ \nabla^2\bm E=\varepsilon_0\mu_0\pdiff{^2\bm E}{t^2} $$となる.最後に,真空中の光速 $c$ は $\varepsilon_0,\mu_0$ を使って $c=1/\sqrt{\varepsilon_0\mu_0}$,すなわち $\varepsilon_0\mu_0=1/c^2$ と表されることを使うと(これは電磁気学における $c$ の定義であり,$\varepsilon_0=8.854\times10^{-12}\,\mathrm{F/m}$,$\mu_0=4\pi\times10^{-7}\,\mathrm{H/m}$ を代入すると実際に $c\approx2.998\times10^8\,\mathrm{m/s}$ になることが確認できる),最終的に
\begin{equation} \nabla^2\bm E=\frac1{c^2}\pdiff{^2\bm E}{t^2} \label{eq:43-wave-E} \end{equation}(証明終わり)
公式43.1 電磁波の波動方程式
真空中では,電場 $\bm E$ は式\eqref{eq:43-wave-E}の波動方程式(wave equation)を満たす.同じ議論を $\bm H$ について行っても(第4式の代わりに第3式から出発すればよい),全く同じ形 $\nabla^2\bm H=\dfrac1{c^2}\pdiff{^2\bm H}{t^2}$ が得られる.つまり,電場も磁場も,真空中を光速 $c$ で伝わる波——電磁波——として振る舞う.マクスウェルが1860年代にこの結果にたどり着き,光そのものが電磁波であると見抜いたことは,物理学史上最大級の発見の1つである.
43.1.2 番外:波動方程式からシュレーディンガー方程式を導く
ここで少し寄り道をしよう.1次元の波動方程式
$$ \pdiff{^2u}{x^2}=\frac1{v^2}\pdiff{^2u}{t^2} $$に,時間について単振動する形の解 $u(x,t)=\psi(x)\cos(\omega t+\theta_0)$($\theta_0$ は初期位相)を代入してみる.$x$ についての偏微分は $\cos(\omega t+\theta_0)$ に影響しないので $\pdiff{^2u}{x^2}=\dfrac{d^2\psi}{dx^2}\cos(\omega t+\theta_0)$ であり,$t$ について2回微分すると $\cos$ を2回微分するので符号が反転して $\pdiff{^2u}{t^2}=-\omega^2\psi\cos(\omega t+\theta_0)$ となる.これらを代入すると,両辺に共通する $\cos(\omega t+\theta_0)$ が割れて,
$$ \frac{d^2\psi}{dx^2}=-\frac{\omega^2}{v^2}\psi $$という,$\psi(x)$ だけについての常微分方程式が得られる.ここで $\omega=2\pi\nu$($\nu$:振動数),$v=\nu\lambda$(波の速さ=振動数×波長,これは高校物理で学ぶ波の基本関係式である)を使うと,
$$ \left(\frac\omega v\right)^2=\left(\frac{2\pi\nu}{\nu\lambda}\right)^2=\frac{4\pi^2}{\lambda^2} $$となる($\nu$ が分子・分母の両方にあるので約分できる).一方,量子力学では,粒子(たとえば電子)は波長 $\lambda$ の波としての性質も持つというド・ブロイ(de Broglie)の関係式 $p=h/\lambda$($h$:プランク定数,$p$:運動量)が成り立つ.粒子の全エネルギー $E$ は運動エネルギーとポテンシャルエネルギー $V$ の和 $E=\dfrac{p^2}{2m}+V$($m$:粒子の質量)で書けるから,これを $p$ について解くと $p=\sqrt{2m(E-V)}$ となり,ド・ブロイの関係式から
$$ \lambda=\frac hp=\frac{h}{\sqrt{2m(E-V)}} $$が得られる.この $\lambda$ を上の $(\omega/v)^2$ の式に代入すると,
$$ \therefore\ \frac{d^2\psi}{dx^2}=-\frac{4\pi^2}{h^2}\cdot2m(E-V)\,\psi $$($4\pi^2/\lambda^2$ に $\lambda^2=h^2/(2m(E-V))$ の逆数 $2m(E-V)/h^2$ を掛けた).両辺に $-\dfrac{h^2}{8\pi^2m}$ を掛けて整理すると(左辺の係数がちょうど $1$ になるように選んだ),
\begin{equation} -\frac{h^2}{8\pi^2m}\pdiff{^2\psi}{x^2}+V\psi=E\psi \label{eq:43-schrodinger-h} \end{equation}という式にたどり着く.
公式43.2 時間に依存しないシュレーディンガー方程式(プランク定数 $h$ による表示)
波動方程式に単振動解を仮定し,ド・ブロイの関係式とエネルギーの表式を組み合わせると,式\eqref{eq:43-schrodinger-h}が得られる.これは時間に依存しないシュレーディンガー方程式(time-independent Schrödinger equation)と呼ばれる,量子力学の基礎方程式そのものである.
なぜ?:この導出は「証明」ではない
ここで行ったのは,シュレーディンガー方程式の厳密な導出ではなく,「粒子が波の性質を持つとしたら,古典的な波動方程式に単振動解を仮定し,ド・ブロイの関係を代入すると,量子力学の基礎方程式と同じ形の式が現れる」というもっともらしさの確認である.実際の量子力学の基礎方程式は,第45章で改めて,演算子(微分演算子としての運動量やエネルギー)という全く別の考え方から導き直される.そこでは,換算プランク定数 $\hbar=h/(2\pi)$ を使った表示 $-\dfrac{\hbar^2}{2m}\pdiff{^2\psi}{x^2}+V\psi=E\psi$ が使われる.$\hbar^2/(2m)=h^2/(8\pi^2m)$($\hbar=h/2\pi$ より $\hbar^2=h^2/4\pi^2$,これを $2m$ で割ると $h^2/(8\pi^2m)$)であることを確かめれば,式\eqref{eq:43-schrodinger-h}が第45章の表示と完全に同じ式であることが分かる.本節の寄り道は,第45章の内容への伏線でもある.
例題43.1 平面波が波動方程式を満たす条件,可視光と電子の波長
(a) $E(x,t)=E_0\cos(kx-\omega t)$($k=2\pi/\lambda$:波数)が1次元の波動方程式 $\pdiff{^2E}{x^2}=\dfrac1{c^2}\pdiff{^2E}{t^2}$ を満たすためには,$\omega,k,c$ の間にどんな関係が必要か.(b) 振動数 $\nu=5.0\times10^{14}\,\mathrm{Hz}$ の可視光の波長を求めよ.(c) 運動エネルギー $100\,\mathrm{eV}$ の電子のド・ブロイ波長を求めよ(電子の質量 $m_e=9.109\times10^{-31}\,\mathrm{kg}$,$1\,\mathrm{eV}=1.602\times10^{-19}\,\mathrm J$,プランク定数 $h=6.626\times10^{-34}\,\mathrm{J\,s}$).
解答 (a) $x$ で2回偏微分すると $\pdiff{^2E}{x^2}=-k^2E_0\cos(kx-\omega t)=-k^2E$,$t$ で2回偏微分すると $\pdiff{^2E}{t^2}=-\omega^2E_0\cos(kx-\omega t)=-\omega^2E$ である.これらを方程式に代入すると
$$ -k^2E=\frac1{c^2}(-\omega^2E)\qquad\therefore\ k^2=\frac{\omega^2}{c^2} $$($E\not\equiv0$ で割った).$k,\omega,c$ はいずれも正の値だから,平方根を取って $\omega=ck$ が条件である.これはまさに図43.1の波が速さ $c$ で伝わるという条件そのものであり,$\omega=ck$(角振動数=伝播速度×波数)は「波の速さ=振動数×波長」という高校物理の関係式 $c=\nu\lambda$ を,$\omega=2\pi\nu$,$k=2\pi/\lambda$ を使って書き直したものに他ならない.
(b) $c=\nu\lambda$ より $\lambda=c/\nu$ である.$c=2.998\times10^8\,\mathrm{m/s}$ を使うと,
$$ \lambda=\frac{2.998\times10^8}{5.0\times10^{14}}=5.996\times10^{-7}\,\mathrm m\approx600\,\mathrm{nm} $$($1\,\mathrm{nm}=10^{-9}\,\mathrm m$).これは橙色(オレンジ色)の可視光の波長域に一致する.
(c) 運動エネルギーは $E_k=\dfrac{p^2}{2m_e}$(ポテンシャル $V=0$ の場合の $E-V$ に相当)だから,$p=\sqrt{2m_eE_k}$ である.まずジュール単位に直すと $E_k=100\times1.602\times10^{-19}=1.602\times10^{-17}\,\mathrm J$.よって
$$ p=\sqrt{2\times9.109\times10^{-31}\times1.602\times10^{-17}}=5.402\times10^{-24}\,\mathrm{kg\,m/s} $$ド・ブロイの関係式 $\lambda=h/p$ より
$$ \lambda=\frac{6.626\times10^{-34}}{5.402\times10^{-24}}=1.227\times10^{-10}\,\mathrm m\approx123\,\mathrm{pm} $$($1\,\mathrm{pm}=10^{-12}\,\mathrm m$).これは原子の大きさ(オングストローム=$100\,\mathrm{pm}$ 程度)と同じ桁であり,だからこそ電子は原子の中で「波」としての性質——定在波としての量子化されたエネルギー準位——を顕著に示すのである(第45章以降で詳しく学ぶ).(sympyとnumpyで検算:(a)は $\omega=ck$ を代入すると恒等的に $0$,(b)(c)は数値計算で上記の値と一致することを確認済み.)
43.2 1次元波動方程式の導出(弦の振動)
前節では電磁気学から波動方程式を導いたが,同じ形の方程式は,ギターやバイオリンの弦の振動からも導くことができる.こちらは高校物理の運動方程式 $F=ma$ だけから出発できるので,波動方程式の意味を最も具体的に理解できる例である.
原点 $0$ と点 $(L,0)$ の2点で固定された弦を考える.弦は一様な線密度 $\rho\,[\mathrm{kg/m}]$(単位長さあたりの質量)と一様な断面積を持ち(ただし,これから立てる1次元近似の運動方程式には断面積は直接には現れない),張力 $T$ は弦のどの位置でも同じ値であるとする.弦を弾くと,弦の各点は上下(図の $u$ 方向)にだけ変位して振動する.時刻 $t$ における位置 $x$ の変位を $u(x,t)$ と書こう.
数学ノート:$\pdiff{^2u}{t^2}$ の書き方について
$u$ は $x$ と $t$ の2変数関数 $u(x,t)$ なので,時刻についての「加速度」は本来 $\dfrac{d^2u}{dt^2}$ ではなく,$x$ を固定した上での時間についての偏微分 $\pdiff{^2u}{t^2}$ と書く.同様に,これから求める合力 $F$ は,弦の各点が動く方向($u$ 軸方向)の力の合力を指す.
弦の微小部分($x$ から $x+\Delta x$ までの区間)を取り出し,そこに働く力を調べよう.この微小部分の質量は $m=\rho\cdot\Delta x$ である.運動方程式 $ma=F$($a$ は $u$ 方向の加速度 $\pdiff{^2u}{t^2}$)から,
\begin{equation} ma=\underbrace{(\rho\cdot\Delta x)}_{m}\cdot\underbrace{\pdiff{^2u}{t^2}}_{a} \label{eq:43-string-newton} \end{equation}である.次に,右辺の合力 $F$ を求める.
導出:合力 $F$ の計算
弦の左端(位置 $x$)では,弦は左下向きに角度 $\theta$ で張力 $T$ を受ける.右端(位置 $x+\Delta x$)では,弦は右上向きに角度 $\theta+\Delta\theta$ で張力 $T$ を受ける.左右の張力の $u$ 方向(鉛直方向)成分は,それぞれ $-T\sin\theta$(左端,下向きなのでマイナス)と $+T\sin(\theta+\Delta\theta)$(右端,上向きなのでプラス)である.したがって,$u$ 方向の合力は
$$ F=T\sin(\theta+\Delta\theta)-T\sin\theta $$である.$\Delta x$ が非常に小さいとき,角度 $\theta,\theta+\Delta\theta$ もまた $0$ に近いので,$\theta\approx0$ のときの近似 $\sin\theta\approx\tan\theta$(これは $\sin\theta=\tan\theta\cos\theta$ で,$\theta\to0$ のとき $\cos\theta\to1$ となることから従う)を使うと,
$$ F\fallingdotseq T\bigl(\tan(\theta+\Delta\theta)-\tan\theta\bigr) $$となる.ここで,$\tan\theta$ は弦の傾き,つまり $u$ を $x$ で微分したもの $\pdiff ux$ に等しい($\tan\theta$ は「$u$ 軸方向にどれだけ登るか/$x$ 軸方向にどれだけ進むか」であり,これはまさに接線の傾きの定義である).位置 $x$ での傾きを $\left(\pdiff ux\right)_x$,位置 $x+\Delta x$ での傾きを $\left(\pdiff ux\right)_{x+\Delta x}$ と書けば,
$$ F=T\left\{\left(\pdiff ux\right)_{x+\Delta x}-\left(\pdiff ux\right)_{x}\right\} $$である.ここで $f(x)=\pdiff ux$(傾きを $x$ の関数とみる)とおくと,上の式は $F=T\{f(x+\Delta x)-f(x)\}$ と書ける.一方,偏微分の定義(極限)を思い出すと,
$$ \pdiff{}{x}f(x)=\lim_{\Delta x\to0}\frac{f(x+\Delta x)-f(x)}{\Delta x} $$であった.これは,$\Delta x$ が十分小さいとき,$f(x+\Delta x)-f(x)\fallingdotseq\Delta x\cdot\pdiff{}{x}f(x)$ という近似が成り立つことを意味する(極限の式の両辺に $\Delta x$ を掛けただけである).$f(x)=\pdiff ux$ なので $\pdiff{}xf(x)=\pdiff{^2u}{x^2}$ であり,したがって
$$ \left(\pdiff ux\right)_{x+\Delta x}-\left(\pdiff ux\right)_{x}=\Delta x\,\pdiff{^2u}{x^2} $$が成り立つ.これを $F$ の式に代入すると,
\begin{equation} F=\Delta x\cdot T\cdot\pdiff{^2u}{x^2} \label{eq:43-string-force} \end{equation}が得られる.
式\eqref{eq:43-string-newton}($ma$ の式)と式\eqref{eq:43-string-force}($F$ の式)は同じ合力を表しているから,両者を等号で結ぶと
$$ \rho\cdot\Delta x\cdot\pdiff{^2u}{t^2}=\Delta x\,T\pdiff{^2u}{x^2} $$両辺を $\rho\cdot\Delta x$ で割ると($\Delta x\ne0$ なので割ってよい),
$$ \pdiff{^2u}{t^2}=\frac T\rho\pdiff{^2u}{x^2} $$ここで $\dfrac T\rho=v^2$($v$ は速さの次元を持つ量として後で確認する)とおくと,
\begin{equation} \pdiff{^2u}{x^2}=\frac1{v^2}\pdiff{^2u}{t^2} \label{eq:43-string-wave} \end{equation}という,43.1節と全く同じ形の波動方程式が得られた.
数学ノート:$v^2=T/\rho$ の次元確認
$T$(張力)の単位は $\mathrm N=\mathrm{kg\cdot m/s^2}$,$\rho$(線密度)の単位は $\mathrm{kg/m}$ である.したがって
$$ [v^2]=\frac{\mathrm{kg\cdot m/s^2}}{\mathrm{kg/m}}=\frac{\mathrm m^2}{\mathrm s^2} $$となり,たしかに速さの2乗の単位($\mathrm{m^2/s^2}$)になっている.方程式 $\pdiff{^2u}{x^2}=\dfrac1{v^2}\pdiff{^2u}{t^2}$ の両辺の次元も確認しておこう.左辺 $\pdiff{^2u}{x^2}$ の次元は $\dfrac{[u]}{\mathrm m^2}$($u$ は変位なので長さの次元),右辺 $\dfrac1{v^2}\pdiff{^2u}{t^2}$ の次元は $\dfrac{\mathrm s^2}{\mathrm m^2}\cdot\dfrac{[u]}{\mathrm s^2}=\dfrac{[u]}{\mathrm m^2}$ となり,両辺の次元は一致する(次元が合わないなら,どこかで計算を間違えている,という検算にも使える).
公式43.3 弦の波動方程式
線密度 $\rho$,張力 $T$ の弦の微小要素に運動方程式 $ma=F$ を適用すると,変位 $u(x,t)$ は波動方程式\eqref{eq:43-string-wave}を満たす.ここで $v=\sqrt{T/\rho}$ は,弦を伝わる波の速さである.
例題43.2 ギターの弦を伝わる波の速さ
あるギターの弦の張力が $T=80\,\mathrm N$,線密度が $\rho=5.0\times10^{-3}\,\mathrm{kg/m}$ であるとき,この弦を伝わる波の速さ $v$ を求めよ.
解答 公式43.3より $v=\sqrt{T/\rho}$ である.数値を代入すると,
$$ v=\sqrt{\frac{80}{5.0\times10^{-3}}}=\sqrt{16000}\approx126.5\,\mathrm{m/s} $$である.弦の長さを $L=0.65\,\mathrm m$(一般的なギターのスケール長)とすれば,基本振動数は 43.4節で学ぶ公式 $\nu_1=\dfrac{v}{2L}$($m=1$ の場合)を先取りして使うと $\nu_1=126.5/(2\times0.65)\approx97\,\mathrm{Hz}$ となり,実際のギターの開放弦(低いE弦,約 $82\,\mathrm{Hz}$)に近い桁の値になる(sympy/numpyで検算:$\sqrt{80/0.005}=126.49\ldots$).
43.3 熱伝導方程式の導出
波動方程式が「振動する」現象を記述するのに対し,本節で導く熱伝導方程式(heat equation)は「拡散する・広がって均される」現象を記述する.細い棒(断面積 $S$)の中の熱の伝わり方を考えよう.棒に沿って座標 $x$ をとり,時刻 $t$,位置 $x$ における温度を $u(x,t)$ とする.熱は温度の高い方から低い方へ,温度分布の勾配(傾き)$\pdiff ux$ に比例して流れる,という経験則(フーリエの熱伝導法則,Fourier's law)を出発点にする.
導出:熱伝導方程式
[i] 流入する熱量 $\Delta Q_1$. 位置 $x+\Delta x$ の断面を通って,時間 $\Delta t$ の間に流入する熱量は,フーリエの熱伝導法則により,その位置での温度勾配 $u_x(x+\Delta x,t)$($u_x=\pdiff ux$ の略記)に比例する.
$$ \Delta Q_1=kS\Delta t\cdot u_x(x+\Delta x,t) $$($k$ は比例定数で熱伝導率と呼ばれる.$S$ は断面積.)
[ii] 流出する熱量 $\Delta Q_2$. 同様に,位置 $x$ の断面を通って流出する熱量は
$$ \Delta Q_2=kS\Delta t\,u_x(x,t) $$である.したがって,時間 $\Delta t$ の間に区間 $[x,x+\Delta x]$ の部分が正味で得た熱量 $\Delta Q$ は,流入から流出を引いたものになる.
$$ \Delta Q=\Delta Q_1-\Delta Q_2=kS\Delta t\bigl\{u_x(x+\Delta x,t)-u_x(x,t)\bigr\} $$ここで,43.2節の弦の議論と全く同じ論法(差分商を偏微分の定義に基づいて近似する)を使うと,
$$ u_x(x+\Delta x,t)-u_x(x,t)=\frac{u_x(x+\Delta x,t)-u_x(x,t)}{\Delta x}\cdot\Delta x\fallingdotseq\pdiff{u_x}x\cdot\Delta x=\Delta x\,\pdiff{^2u}{x^2} $$(最後の等号は $\pdiff{}x u_x=\pdiff{}x\pdiff ux=\pdiff{^2u}{x^2}$ による).したがって,
\begin{equation} \Delta Q\fallingdotseq kS\Delta t\cdot\Delta x\,\pdiff{^2u}{x^2} \label{eq:43-heat-Q1} \end{equation}となる.
[iii] 熱量と温度上昇の関係. 一方,この熱量 $\Delta Q$ によって区間 $[x,x+\Delta x]$ の部分の温度が $\Delta u$ だけ上昇したとすると,物理でおなじみの関係式「熱量=質量×比熱×温度上昇」($Q=mC\Delta u$,$C$ は比熱)が成り立つ.この区間の質量は,体積密度を $\rho$ として $m=\rho S\Delta x$(体積 $S\Delta x$ に密度を掛けたもの)である.また,温度上昇 $\Delta u$ は,時間 $\Delta t$ の間の変化なので,$\Delta u\fallingdotseq\Delta t\,\pdiff ut$ と近似できる(これも偏微分の定義そのものである).よって
\begin{equation} \Delta Q=m\cdot C\cdot\Delta u=\rho S\Delta x\cdot C\cdot\Delta u\fallingdotseq\rho SC\Delta x\cdot\Delta t\,\pdiff ut \label{eq:43-heat-Q2} \end{equation}である.
[iv] 2通りに表した $\Delta Q$ を等号で結ぶ. 式\eqref{eq:43-heat-Q1}と式\eqref{eq:43-heat-Q2}はどちらも同じ熱量 $\Delta Q$ を表しているから,
$$ kS\Delta t\,\Delta x\pdiff{^2u}{x^2}=\rho SC\Delta x\,\Delta t\pdiff ut $$両辺を $S\Delta x\Delta t$ で割ると(共通因子なので消える),
$$ k\pdiff{^2u}{x^2}=\rho C\pdiff ut $$最後に $\dfrac{k}{\rho C}=a$($a$ は温度伝導率と呼ばれる)とおいて整理すると,
\begin{equation} \pdiff ut=a\pdiff{^2u}{x^2} \label{eq:43-heat-eq} \end{equation}が導かれる.
(証明終わり)
公式43.4 熱伝導方程式
断面積 $S$,密度 $\rho$,比熱 $C$,熱伝導率 $k$ の棒の微小区間で熱収支を考えると,温度分布 $u(x,t)$ は熱伝導方程式\eqref{eq:43-heat-eq}を満たす.ここで $a=k/(\rho C)$ は温度伝導率である.波動方程式\eqref{eq:43-string-wave}が時間について2階の偏微分方程式だったのに対し,熱伝導方程式は時間について1階の偏微分方程式である点に注意しよう.この違いが,43.4節以降で見るように,解の性質(振動し続けるか,指数的に減衰するか)を決定的に分ける.
例題43.3 熱伝導方程式の定常解
棒の両端の温度が時間によらず一定値 $u(0,t)=u_1$,$u(L,t)=u_2$ に保たれているとき,十分に時間が経った後の定常状態(温度分布がもはや時間変化しない状態,$\pdiff ut=0$)における温度分布 $u(x)$ を求めよ.具体的に $u_1=20\,{}^\circ\mathrm C$,$u_2=80\,{}^\circ\mathrm C$,$L=2\,\mathrm m$ のとき,中点 $x=1\,\mathrm m$ での温度を求めよ.
解答 定常状態では $\pdiff ut=0$ だから,熱伝導方程式\eqref{eq:43-heat-eq}は $a\pdiff{^2u}{x^2}=0$ となる.$a\ne0$ なので,これは常微分方程式 $\dfrac{d^2u}{dx^2}=0$($t$ に依存しなくなったので $u=u(x)$ の常微分方程式になる)と同じであり,2回積分すると一般解は
$$ u(x)=Ax+B $$($A,B$ は積分定数)という1次関数(線形分布)になる.境界条件 $u(0)=u_1$ より $B=u_1$,$u(L)=u_2$ より $AL+u_1=u_2$,すなわち $A=(u_2-u_1)/L$ である.したがって
$$ u(x)=u_1+\frac{u_2-u_1}{L}x $$数値を代入すると,$u(x)=20+\dfrac{80-20}{2}x=20+30x$ となる.$x=1$ での温度は
$$ u(1)=20+30\times1=50\,{}^\circ\mathrm C $$である(中点なので,直感どおり両端の平均 $(20+80)/2=50\,{}^\circ\mathrm C$ に一致する——これは分布が線形だから成り立つ性質である).熱伝導方程式の定常解が常に線形分布(グラフでは直線)になることは,「熱は温度勾配に比例して流れる」という物理法則の帰結として自然である:定常状態では,どの断面を通る熱の流れも一定でなければならず(そうでないとどこかに熱がたまって温度が変化してしまう),流れが勾配に比例するのだから,勾配自体が一定,つまりグラフが直線になるのである.(sympyで検算:$u(x)=20+30x$ が $u''=0$ を満たし,$u(0)=20,u(2)=80,u(1)=50$ となることを確認済み.)
43.4 1次元波動方程式の変数分離解法
43.1節・43.2節で波動方程式そのものを導いたので,いよいよ具体的な境界値問題を解く.ここで使う道具が変数分離法(method of separation of variables)である.
定義43.1 変数分離法
2変数関数 $u(x,t)$ についての偏微分方程式を解くとき,解を $x$ だけの関数 $X(x)$ と $t$ だけの関数 $T(t)$ の積 $u(x,t)=X(x)T(t)$ という特別な形に限定して探す方法を,変数分離法という.この形を偏微分方程式に代入すると,多くの場合,方程式が「$x$ だけを含む部分」と「$t$ だけを含む部分」に分離できる.両辺は独立変数 $x,t$ をそれぞれ含むので,両辺が定数(分離定数)に等しくなるしかない.こうして,もとの偏微分方程式が,$X(x)$ についての常微分方程式と $T(t)$ についての常微分方程式という,2つの(より簡単な)常微分方程式に分解される.
43.4.1 問題設定と変数分離
次の境界値問題を考えよう.長さ $L$ の弦(両端 $x=0,L$ で固定)が,時刻 $t=0$ に図43.5のような三角形の形に変位させられ,初速度 $0$ の状態から手を離されたとする.
例43.1 1次元波動方程式:三角形分布の初期変位(ノートの例題)
関数 $u(x,t)$ について,波動方程式
$$ \pdiff{^2u}{t^2}=v^2\pdiff{^2u}{x^2}\qquad(0\lt x\lt L,\ 0\lt t) $$を,初期条件
$$ u(x,0)=\begin{cases}\dfrac{3h}Lx&\left(0\le x\le\dfrac L3\right)\\[6pt]\dfrac{3h}{2L}(L-x)&\left(\dfrac L3\lt x\le L\right)\end{cases}\, ,\qquad \pdiff ut(x,0)=0 $$境界条件 $u(0,t)=u(L,t)=0$ のもとで解け.
以下,この問題を解いていく.
解答(1)変数分離 変数分離法により,$u(x,t)=X(x)T(t)$ とおく(本来は $u=X(x)T(t)$ の $T$ と,張力の $T$ は別の文字だが,ノートの表記に合わせてこのまま用いる.文脈で混同しないよう注意する).これを波動方程式に代入すると,$\pdiff{^2u}{t^2}=X\ddot T$,$\pdiff{^2u}{x^2}=X''T$($\ddot T=d^2T/dt^2$,$X''=d^2X/dx^2$)であるから,
$$ X\ddot T=v^2X''T $$両辺を $v^2XT$ で割ると($X,T$ がともに恒等的に $0$ でないと仮定する),
$$ \frac{X''}X=\frac{\ddot T}{v^2T} $$という形になる.左辺 $X''/X$ は $x$ だけの関数($t$ を含まない)であり,右辺 $\ddot T/(v^2T)$ は $t$ だけの関数($x$ を含まない)である.なぜこの式の両辺が定数でなければならないのか,ここで確認しておこう.左辺は $x$ だけの関数なので,$t$ をどう動かしても左辺の値は変わらない.ところが両辺は等しいのだから,$t$ を動かしたときの右辺の値も変わってはならない——つまり右辺は,実は $t$ によらない($t$ を含む式なのに,$t$ を動かしても値が変わらないのだから,結局 $x$ だけの何かに等しいはずである).同様に,今度は $x$ を動かしてみると,右辺($t$ だけの関数のはず)はやはり変わらないのだから,左辺もまた $x$ によらない.結局,この共通の値は,$x,t$ のどちらを動かしても変化しない——すなわち定数でなければならない.この定数を分離定数(separation constant)と呼び,$C_0$ と書くことにする.
$$ \frac{X''}X=\frac{\ddot T}{v^2T}=C_0\qquad(C_0:\text{分離定数}) $$43.4.2 固有値問題の場合分け
数学ノート:なぜ非自明解の条件を調べるのか
$X''=C_0X$ の一般解は,$C_0$ の符号に応じて指数関数($C_0\gt0$),1次関数($C_0=0$),三角関数($C_0\lt0$)のいずれかになる.それぞれの場合で境界条件 $X(0)=X(L)=0$ を課すと,未定係数についての斉次連立1次方程式(右辺がすべて $0$ の連立方程式)が出てくる.第8章・第9章で学んだように,斉次連立1次方程式が「係数がすべて $0$」という自明な解以外の解(非自明解)を持つのは,係数行列の行列式が $0$ になるとき,かつそのときに限る.ここでは行列式を陽に書き下さずに直接計算するが,各場合で「非自明解が存在するか」を調べるという作業は,まさにこの行列式の条件を確認していることに他ならない.なぜ非自明解にこだわるかというと,$X(x)\equiv0$(自明解)では $u=XT\equiv0$ となってしまい,何も振動していない「解」しか得られないからである.
[1] $X''=C_0X$ について,非自明な $X$ が存在する $C_0$ を探す.
(i) $C_0\gt0$ のとき.特性方程式(第35章)の解は $\pm\sqrt{C_0}$ なので,一般解は
$$ X(x)=A_1e^{\sqrt{C_0}x}+B_1e^{-\sqrt{C_0}x} $$境界条件 $X(0)=0$ より $A_1+B_1=0$.境界条件 $X(L)=0$ より $A_1e^{\sqrt{C_0}L}+B_1e^{-\sqrt{C_0}L}=0$.この2つの条件を同時に満たす $A_1,B_1$ を求めよう.$B_1=-A_1$ を2つ目の式に代入すると $A_1(e^{\sqrt{C_0}L}-e^{-\sqrt{C_0}L})=0$ となるが,$L\gt0$,$\sqrt{C_0}\gt0$ より $e^{\sqrt{C_0}L}\ne e^{-\sqrt{C_0}L}$(指数関数は狭義単調増加だから,異なる指数には異なる値が対応する)なので,かっこの中は $0$ にならない.したがって $A_1=0$,ゆえに $B_1=0$ しかない.つまりこの場合をみたすのは自明解 $X(x)\equiv0$ だけであり,不適である.
(ii) $C_0=0$ のとき.$X''=0$ より,一般解は
$$ X(x)=A_2x+B_2 $$境界条件 $X(0)=0$ より $B_2=0$.境界条件 $X(L)=0$ より $A_2L+B_2=0$,$B_2=0$ を代入すると $A_2L=0$,$L\ne0$ より $A_2=0$.やはり自明解しかなく,不適である.
(iii) $C_0\lt0$ のとき.$C_0=-\omega^2$($\omega\gt0$)とおくと,特性方程式の解は $\pm i\omega$ なので,一般解は三角関数になる(第35章).
$$ X(x)=A_3\cos\omega x+B_3\sin\omega x $$境界条件 $X(0)=0$ より $A_3=0$.境界条件 $X(L)=0$ より $A_3\cos\omega L+B_3\sin\omega L=0$,$A_3=0$ を代入すると $B_3\sin\omega L=0$ となる.ここで,もし $B_3=0$ ならまたも自明解になってしまうので,非自明解を得るには $B_3\ne0$ でなければならず,そのためには
$$ \sin\omega L=0 $$が必要である.$\omega\gt0$,$L\gt0$ なので $\omega L\gt0$,したがって $\omega L=m\pi$($m=1,2,3,\ldots$;$m=0$ は $\omega=0$ となり $C_0=0$ に戻ってしまうので除く)である.こうして,非自明解が存在するのは $C_0=-\omega^2=-(m\pi/L)^2$($m=1,2,\ldots$)というとびとびの値のときだけであることが分かった.このように,境界条件によって特定の値だけが許される定数 $C_0$(あるいは $\omega$)を固有値(eigenvalue),対応する $X(x)$ を固有関数(eigenfunction)と呼ぶ.固有値・固有関数という言葉は,第12章で学んだ行列の固有値・固有ベクトルと同じ発想であり,「特別な入力に対してだけ,出力が入力の定数倍になる」という構造が共通している.
このとき,固有関数は
$$ X(x)=B_3\sin\frac{m\pi}Lx\qquad(m=1,2,3,\ldots) $$である.
[2] $T(t)$ を求める. 固有値 $C_0=-(m\pi/L)^2$ が定まったので,これを $T$ の方程式に代入する.
$$ \ddot T=C_0v^2T=-\left(\frac{m\pi}Lv\right)^2T $$これは単振動の方程式(第35章)だから,一般解は
$$ T(t)=C_1\cos\frac{m\pi}Lvt+D_1\sin\frac{m\pi}Lvt $$ここで初期条件 $\pdiff ut(x,0)=0$(初速度が $0$)を使う.$\pdiff ut=X(x)\dot T(t)$ で,$\dot T(t)=-C_1\dfrac{m\pi}Lv\sin\dfrac{m\pi}Lvt+D_1\dfrac{m\pi}Lv\cos\dfrac{m\pi}Lvt$ だから,$\dot T(0)=D_1\dfrac{m\pi}Lv=0$,したがって($m,\pi,L,v$ はすべて $0$ でないので)$D_1=0$.よって
$$ T(t)=C_1\cos\frac{m\pi}Lvt $$43.4.3 重ね合わせとフーリエ係数の決定
[1],[2]から,各 $m=1,2,3,\ldots$ に対して,波動方程式と境界条件を満たす解
$$ U_m(x,t)=B_3C_1\sin\frac{m\pi}Lx\cdot\cos\frac{m\pi}Lvt $$が得られる($B_3C_1$ をまとめて $D_m$ とおく:正弦のフーリエ級数を扱うことを見越した記号である).しかし,$U_m(x,t)$ はただ1つの $m$ に対応する解であり,一般には三角形分布の初期条件 $u(x,0)$ を満たさない.ここで重要な性質を使う.
定理43.1 重ね合わせの原理
波動方程式 $\pdiff{^2u}{t^2}=v^2\pdiff{^2u}{x^2}$ は線形偏微分方程式である($u$ について1次式であり,$u$ の積や $u$ を含む関数などが現れない).したがって,$U_1,U_2,\ldots$ がそれぞれこの方程式の解であれば,それらの無限1次結合(線形結合ともいう)
$$ u=\sum_{m=1}^\infty D_mU_m $$もまた解になる.なぜなら,方程式の左辺・右辺はどちらも $u$ について1次の微分演算であり,和の微分は微分の和に等しい(線形性)から,各項がそれぞれ方程式を満たしていれば,その和も方程式を満たすからである.
この重ね合わせの原理により,一般解は
$$ u(x,t)=\sum_{m=1}^\infty U_m(x,t)=\sum_{m=1}^\infty D_m\sin\frac{m\pi}Lx\cos\frac{m\pi}Lvt $$という無限級数の形になる.残る仕事は,初期条件 $u(x,0)=(\text{三角形分布})$ を満たすように,係数 $D_m$ を決定することである.$t=0$ を代入すると $\cos0=1$ なので,
$$ u(x,0)=\sum_{m=1}^\infty D_m\sin\frac{m\pi}Lx $$となる.これはまさに,$u(x,0)$ を区間 $[0,L]$ 上のフーリエ正弦級数(第28章・第29章)に展開した式である.フーリエ正弦級数の係数公式(正弦関数の直交性 $\displaystyle\int_0^L\sin\frac{m\pi}Lx\sin\frac{m'\pi}Lx\,dx=\dfrac L2\delta_{mm'}$ から導かれる)
$$ D_m=\frac2L\int_0^Lu(x,0)\sin\frac{m\pi}Lx\,dx $$を使う.初期条件が2つの式に分かれているので,積分区間も分けて計算する.
$$ D_m=\frac2L\left\{\int_0^{L/3}\frac{3h}Lx\sin\frac{m\pi}Lx\,dx+\int_{L/3}^L\frac{3h}{2L}(L-x)\sin\frac{m\pi}Lx\,dx\right\} $$それぞれの積分は,部分積分(第5章)を2回繰り返す標準的な計算である($\displaystyle\int x\sin ax\,dx=-\dfrac xa\cos ax+\dfrac1{a^2}\sin ax$ の形の公式を使う).
$$ \frac{3h}L\int x\sin\frac{m\pi}Lx\,dx=\frac{3h}L\left(-\frac L{m\pi}x\cos\frac{m\pi}Lx+\frac{L^2}{m^2\pi^2}\sin\frac{m\pi}Lx\right) $$ $$ \int\frac{3h}{2L}(L-x)\sin\frac{m\pi}Lx\,dx=-\frac{3hL}{2m\pi}\cos\frac{m\pi}Lx-\frac{3h}{2L}\left(-\frac L{m\pi}x\cos\frac{m\pi}Lx+\frac{L^2}{m^2\pi^2}\sin\frac{m\pi}Lx\right) $$これらの原始関数に積分区間の上端・下端を代入し(三角関数の値 $\sin(m\pi/3),\cos(m\pi/3),\sin m\pi=0,\cos m\pi=(-1)^m$ などを使って整理すると,途中の多くの項が打ち消し合う),最終的に
\begin{equation} D_m=\frac{9h}{m^2\pi^2}\sin\frac{m\pi}3 \label{eq:43-Dm} \end{equation}という簡潔な形にまとまる.したがって,求める解は
\begin{equation} u(x,t)=\sum_{m=1}^\infty\frac{9h}{m^2\pi^2}\sin\frac{m\pi}3\cdot\sin\frac{m\pi}Lx\cos\frac{m\pi}Lvt \label{eq:43-wave-triangle-sol} \end{equation}(例43.1の解答終わり)
数学ノート:$\sin(m\pi/3)$ の値と,消える項
$\sin(m\pi/3)$ は $m$ を $6$ で割った余りによって周期的に符号や値が変わる:$m=1$ で $\sin(\pi/3)=\sqrt3/2$,$m=2$ で $\sin(2\pi/3)=\sqrt3/2$,$m=3$ で $\sin\pi=0$,$m=4$ で $\sin(4\pi/3)=-\sqrt3/2$,$m=5$ で $\sin(5\pi/3)=-\sqrt3/2$,$m=6$ で $\sin(2\pi)=0$,そして $m=7$ からまた同じ周期を繰り返す.つまり,$3$ の倍数の $m$ に対応する項($D_3,D_6,D_9,\ldots$)はちょうど $0$ になり,級数から消える.これは初期条件の三角形の「山」の頂点がちょうど $x=L/3$(区間の $1/3$ の位置)にあることと関係している(sympyでの検算は本節末の演習と例題43.4で行う).
図43.6に,この解を実際に数値計算して,弦の形の時間変化を追った様子を示す($L=1,h=1,v=1$ として無次元化し,級数を $m=1$〜$400$ まで足し合わせた).基本振動の周期は $T_{\mathrm{period}}=2L/v$ である.
例題43.4 初期条件が単振動1つで済む場合
長さ $L$,両端固定の弦が,初期条件 $u(x,0)=h\sin\dfrac{\pi x}L$,$\pdiff ut(x,0)=0$ から出発するとき,波動方程式の解を求めよ.
解答 43.4.3項までの議論はそのまま使えて,一般解は $u(x,t)=\displaystyle\sum_{m=1}^\infty D_m\sin\dfrac{m\pi}Lx\cos\dfrac{m\pi}Lvt$($D_m$ はフーリエ正弦級数の係数)である.今回の初期条件は,最初から $u(x,0)=h\sin(\pi x/L)$ という1つの固有関数($m=1$ の場合)そのものになっている.正弦関数の直交性 $\displaystyle\int_0^L\sin\frac{\pi x}L\sin\frac{m\pi x}Ldx=\frac L2\delta_{1m}$($m=1$ のときだけ $L/2$,それ以外は $0$)を係数公式に使うと,
$$ D_m=\frac2L\int_0^Lh\sin\frac{\pi x}L\sin\frac{m\pi x}Ldx=\begin{cases}h&(m=1)\\0&(m\ne1)\end{cases} $$となる.つまり無限級数のうち $m=1$ の項だけが生き残り,求める解は
$$ u(x,t)=h\sin\frac{\pi x}L\cos\frac{\pi vt}L $$という,たった1項の式になる.これは,最初から「形」が固有関数と一致している初期条件を与えると,時間発展してもその「形」($\sin(\pi x/L)$)を保ったまま,振幅だけが $\cos(\pi vt/L)$ に従って振動する,という美しい結果を示している.一般の初期条件(例43.1の三角形分布のような)は,このような単純な固有振動の無限個の重ね合わせだと考えることができる.(sympyで検算:$m=1$ のとき $D_1=h$,$m\ge2$ のとき $D_m=0$ となることを直接積分して確認済み.)
43.5 2次元波動方程式の変数分離解法
前節では弦(1次元)を扱ったが,同じ手法は,太鼓の膜のような2次元の振動にもそのまま拡張できる.変数分離法の考え方は「独立変数の数だけ関数を掛け合わせる」というだけなので,2変数が3変数($x,y,t$)に増えても本質的には同じ手順を繰り返すだけである.
例43.2 2次元波動方程式:山型の初期変位(ノートの例題)
関数 $u(x,y,t)$ について,波動方程式
$$ \pdiff{^2u}{t^2}=v^2\left(\pdiff{^2u}{x^2}+\pdiff{^2u}{y^2}\right)\qquad(0\lt x\lt1,\ 0\lt y\lt1,\ 0\lt t) $$を,初期条件 $u(x,y,0)=x(1-x)y(1-y)$,$\pdiff ut(x,y,0)=0$,境界条件 $u(0,y,t)=u(1,y,t)=u(x,0,t)=u(x,1,t)=0$(正方形の膜の縁が固定されている)のもとで解け.
解答(1)変数分離 $u(x,y,t)=X(x)Y(y)T(t)$ とおく.波動方程式に代入すると,$\pdiff{^2u}{t^2}=XY\ddot T$,$\pdiff{^2u}{x^2}=X''YT$,$\pdiff{^2u}{y^2}=XY''T$ であるから,
$$ XY\ddot T=v^2(X''YT+XY''T) $$両辺を $v^2XYT$ で割ると,
$$ \frac{\ddot T}{v^2T}=\frac{X''}X+\frac{Y''}Y $$左辺は $t$ だけの関数,右辺は $x,y$ だけの関数だから,この式が恒等的に成り立つには,両辺がある定数($-\omega^2$ とおく)に等しくなければならない(なぜそう言えるかは43.4.1項で確認した理由と全く同じである:一方を固定してもう一方を動かしても値が変わらないことから,両辺は共通の定数でなければならない).
$$ \frac{\ddot T}{v^2T}=\frac{X''}X+\frac{Y''}Y=-\omega^2 $$さらに,右辺 $\dfrac{X''}X+\dfrac{Y''}Y$ 自身も,$\dfrac{X''}X$ が $x$ だけの関数,$\dfrac{Y''}Y$ が $y$ だけの関数の和だから,同じ理屈で各項がそれぞれ定数でなければならない.そこで $\dfrac{X''}X=-\omega_1^2$,$\dfrac{Y''}Y=-\omega_2^2$($\omega_1^2+\omega_2^2=\omega^2$)とおくと,3つの常微分方程式
$$ X''=-\omega_1^2X,\qquad Y''=-\omega_2^2Y,\qquad \ddot T=-v^2\omega^2T $$に分解される(前節の1次元の場合よりも分離定数が1つ増えただけで,考え方は完全に同じである).
解答(2)固有値問題 境界条件 $X(0)=X(1)=0$ のもとで $X''=-\omega_1^2X$ を解く.一般解は $X(x)=A_1\cos\omega_1x+B_1\sin\omega_1x$.$X(0)=A_1=0$.$X(1)=B_1\sin\omega_1=0$ で,非自明解 $B_1\ne0$ を得るには $\sin\omega_1=0$,すなわち $\omega_1=m\pi$($m=1,2,\ldots$)でなければならない(前節43.4.2項の議論と全く同じ理由による).よって
$$ X(x)=B_1\sin m\pi x\qquad(m=1,2,\ldots) $$全く同様に,境界条件 $Y(0)=Y(1)=0$ のもとで $Y''=-\omega_2^2Y$ を解くと,
$$ Y(y)=B_2\sin n\pi y\qquad(n=1,2,\ldots) $$($\omega_2=n\pi$)が得られる.最後に $T$ の方程式 $\ddot T=-v^2\omega^2T$(単振動)を解く.初期条件 $\pdiff ut(x,y,0)=0$ から,前節と同様の議論で $\sin$ の係数が $0$ になり,
$$ T(t)=A_3\cos v\omega t $$となる.ここで $\omega=\sqrt{\omega_1^2+\omega_2^2}=\sqrt{m^2\pi^2+n^2\pi^2}=\sqrt{m^2+n^2}\,\pi$ である.
解答(3)重ね合わせとフーリエ係数 以上をまとめると,各 $m,n=1,2,\ldots$ に対する解は
$$ U_{mn}(x,y,t)=B_1B_2A_3\sin m\pi x\sin n\pi y\cos v\omega t $$である($B_1B_2A_3=D_{mn}$ とおく).定理43.1(重ね合わせの原理)は,$x,y,t$ の3変数になっても全く同じ理由で成り立つので,一般解は無限二重級数
$$ u(x,y,t)=\sum_{n=1}^\infty\sum_{m=1}^\infty D_{mn}\sin m\pi x\sin n\pi y\cos v\sqrt{m^2+n^2}\,\pi t $$になる.$t=0$ を代入すると,初期条件から
$$ u(x,y,0)=\sum_{n=1}^\infty\sum_{m=1}^\infty D_{mn}\sin m\pi x\sin n\pi y=x(1-x)y(1-y) $$これは2重フーリエ正弦級数(第28章・第29章の1変数の場合を $x,y$ それぞれの方向に独立に適用したもの)であり,係数は次のように,$x,y$ それぞれについての1変数の係数公式を掛け合わせた形になる.
$$ D_{mn}=\frac4{L_1L_2}\int_0^{L_2}\int_0^{L_1}u(x,y,0)\sin\frac{m\pi}{L_1}x\sin\frac{n\pi}{L_2}y\,dx\,dy $$今回は $L_1=L_2=1$ で,$u(x,y,0)=x(1-x)y(1-y)$ が $x$ の部分と $y$ の部分の積の形をしているので,2重積分は2つの1重積分の積に分解できる.
$$ D_{mn}=4\int_0^1x(1-x)\sin m\pi x\,dx\int_0^1y(1-y)\sin n\pi y\,dy $$この形の積分は,43.4節と同じく部分積分を2回繰り返せば計算できる.
$$ \int_0^1(x-x^2)\sin m\pi x\,dx=\left[(x^2-x)\frac{\cos m\pi x}{m\pi}+\frac{\sin m\pi x}{m\pi}-\frac2{m^2\pi^2}\left(x\sin m\pi x+\frac{\cos m\pi x}{m\pi}\right)\right]_0^1 $$両端点 $x=0,1$ を代入すると($x^2-x$ は $x=0,1$ でどちらも $0$,$\sin m\pi=0$ であることに注意),最終的に
$$ \int_0^1(x-x^2)\sin m\pi x\,dx=-\frac2{m^2\pi^2}\left(\frac{(-1)^m}{m\pi}-\frac1{m\pi}\right)=\frac{2\bigl(1-(-1)^m\bigr)}{m^3\pi^3} $$となる($y$ 方向の積分も文字を $n$ に変えるだけで全く同じ形になる).したがって
$$ D_{mn}=4\cdot\frac{2\bigl(1-(-1)^m\bigr)}{m^3\pi^3}\cdot\frac{2\bigl(1-(-1)^n\bigr)}{n^3\pi^3}=\frac{16}{\pi^6}\cdot\frac{1-(-1)^m}{m^3}\cdot\frac{1-(-1)^n}{n^3} $$が得られ,求める解は
\begin{equation} u(x,y,t)=\frac{16}{\pi^6}\sum_{n=1}^\infty\sum_{m=1}^\infty\frac{1-(-1)^m}{m^3}\cdot\frac{1-(-1)^n}{n^3}\sin m\pi x\sin n\pi y\cos v\sqrt{m^2+n^2}\,\pi t \label{eq:43-wave2d-sol} \end{equation}である(例43.2の解答終わり).
注意:ノートの誤記について
この最終解(式\eqref{eq:43-wave2d-sol})を導く直前の $D_{mn}$(あるいは $b_{mn}$)の式では分母が $m^3,n^3$ であったが,ノートの最終行では分母が $m,n$(3乗が付いていない)と書かれている.sympyで独立に積分を計算し直すと,$D_{mn}$ の式(分母 $m^3,n^3$)が正しく,最終行の分母 $m,n$ は書き写しの際の誤記であることが確認できる.本章では正しい形(分母 $m^3,n^3$)で統一した.
図43.8に,この解を数値計算した様子を示す(ノートと同じく,級数を $m,n=1$〜$15$ で打ち切ってシミュレーションしている).
例題43.5 初期条件が単振動1つで済む場合(2次元)
単位正方形 $0\lt x\lt1,\ 0\lt y\lt1$ の膜が,初期条件 $u(x,y,0)=\sin\pi x\sin\pi y$,$\pdiff ut(x,y,0)=0$ から出発するとき,波動方程式の解を求めよ.
解答 例43.2までの結果から,一般解は $u(x,y,t)=\displaystyle\sum_{n,m}D_{mn}\sin m\pi x\sin n\pi y\cos v\sqrt{m^2+n^2}\,\pi t$ である.今回の初期条件 $\sin\pi x\sin\pi y$ は,$m=n=1$ の固有関数そのものなので,2重フーリエ正弦級数の直交性($x$ 方向・$y$ 方向それぞれで例題43.4と同じ直交性が働く)から,$D_{11}=1$,それ以外の $D_{mn}=0$ となる.したがって解は
$$ u(x,y,t)=\sin\pi x\sin\pi y\cos\bigl(v\sqrt2\,\pi t\bigr) $$である($\omega=\sqrt{1^2+1^2}\,\pi=\sqrt2\,\pi$).例題43.4(1次元)と同様,初期条件が固有関数そのものであれば,時間発展してもその「形」を保ったまま振動する.(sympyで検算:$m=n=1$ で $b_{11}=1$,他の $m,n\,(\le3)$ ではすべて $0$ となることを直接積分して確認済み.)
43.6 1次元熱伝導方程式の変数分離解法とシミュレーション
波動方程式と全く同じ手順(変数分離 → 固有値問題 → フーリエ係数の決定)を,今度は熱伝導方程式に適用する.方程式の形は違うが,$X(x)$ についての固有値問題は波動方程式の場合と全く同じ形になる——境界条件が同じ「両端 $u=0$」だからである.違いが現れるのは $T(t)$ の側であり,そこに熱伝導方程式の本質的な性質(振動せず,指数的に減衰する)が表れる.
例43.3 1次元熱伝導方程式:矩形パルスの初期温度分布(ノートの例題)
関数 $U(x,t)$ について,熱伝導方程式
$$ \pdiff ut=a\pdiff{^2u}{x^2}\qquad(0\lt x\lt L,\ 0\lt t)\qquad(a\gt0,\ T_0\gt0) $$を,初期条件
$$ U(x,0)=\begin{cases}T_0&\left(\dfrac L2\le x\le\dfrac34L\right)\\0&\left(0\le x\lt\dfrac L2,\ \dfrac34L\lt x\le L\right)\end{cases} $$境界条件 $U(0,t)=U(L,t)=0$ のもとで解け.
解答(1)変数分離と固有値問題 $U(x,t)=X(x)T(t)$ とおく.熱伝導方程式に代入すると $X\dot T=aX''T$,両辺を $aXT$ で割ると
$$ \frac{X''}X=\frac{\dot T}{aT}=C_0\qquad(C_0:\text{定数}) $$境界条件 $X(0)=X(L)=0$ のもとで非自明解が存在するのは,43.4.2項と全く同じ議論により $C_0\lt0$ のときだけである($C_0\ge0$ だと自明解しか出ない).$C_0=-\omega^2$($\omega\gt0$)とおくと,
$$ \frac{X''}X=\frac{\dot T}{aT}=-\omega^2 $$[1] $X''=-\omega^2X$ について,一般解は $X(x)=A_1\cos\omega x+B_1\sin\omega x$.境界条件から $X(0)=A_1=0$ かつ $X(L)=B_1\sin\omega L=0$.非自明解を得るには $\sin\omega L=0$,すなわち $\omega L=n\pi$($n\in\N$)である.
[2] $\dot T=-\omega^2aT$ について.これは1階の常微分方程式であり,変数分離形(第33章)としてただちに解ける.$\dfrac{dT}T=-\omega^2a\,dt$ を積分すると $\ln|T|=-\omega^2at+(\text{定数})$,したがって
$$ T(t)=A_2e^{-\omega^2at}=A_2e^{-\frac{n^2\pi^2}{L^2}at} $$という指数関数的に減衰する解になる.これが,波動方程式の $T(t)=C_1\cos(\cdots)$(振動し続ける)との決定的な違いである——熱伝導方程式は時間について1階なので,$T$ の方程式は単振動ではなく指数関数の減衰になるのである.
解答(2)重ね合わせとフーリエ係数 [1],[2]から,各 $n=1,2,\ldots$ に対する解は
$$ U_n(x,t)=B_1A_2\,e^{-\frac{n^2\pi^2}{L^2}at}\sin\frac{n\pi}Lx $$である($B_1A_2=b_n$ とおく).熱伝導方程式も線形なので,重ね合わせの原理(定理43.1)が同様に成り立ち,一般解は
$$ U(x,t)=\sum_{n=1}^\infty U_n(x,t)=\sum_{n=1}^\infty b_n\,e^{-\frac{n^2\pi^2}{L^2}at}\sin\frac{n\pi}Lx $$となる.$t=0$ を代入すると $U(x,0)=\displaystyle\sum_n b_n\sin\dfrac{n\pi}Lx$ となり,これはフーリエ正弦級数の係数公式
$$ b_n=\frac2L\int_0^LU(x,0)\sin\frac{n\pi}Lx\,dx $$で決まる.初期条件は $[L/2,3L/4]$ の外側で $0$ なので,積分区間はそこだけに絞られる.
$$ b_n=\frac2L\int_{L/2}^{3L/4}T_0\sin\frac{n\pi}Lx\,dx=\left[-\frac{2T_0}L\cdot\frac L{n\pi}\cos\frac{n\pi}Lx\right]_{L/2}^{3L/4}=-\frac{2T_0}{n\pi}\left(\cos\frac{3n}4\pi-\cos\frac n2\pi\right) $$したがって,求める解は
\begin{equation} U(x,t)=\sum_{n=1}^\infty\frac{2T_0}{n\pi}\left(\cos\frac{n\pi}2-\cos\frac{3n}4\pi\right)e^{-\frac{n^2\pi^2}{L^2}at}\sin\frac{n\pi}Lx \label{eq:43-heat1d-sol} \end{equation}である(例43.3の解答終わり).
公式43.5 波動方程式と熱伝導方程式の解の質的な違い
同じ「両端固定」の境界条件のもとで変数分離すると,空間方向の固有関数 $\sin(n\pi x/L)$ と固有値 $\omega_n=n\pi/L$ は,波動方程式でも熱伝導方程式でも共通である.違いは時間発展の部分だけに現れる:波動方程式では $T_n(t)=\cos(n\pi vt/L)$(永久に振動し続ける),熱伝導方程式では $T_n(t)=e^{-n^2\pi^2at/L^2}$($n$ が大きいほど速く減衰する).この違いは,方程式の時間微分の階数(2階か1階か)の違いに由来する.物理的にも自然である:弦はエネルギーを失わなければ振動し続けるが,熱は温度差がある限り高温部から低温部へ流れ続け,最終的に一様な温度(この境界条件のもとでは $U\equiv0$)に近づいていく.
図43.10に,実際の時間発展を数値計算した様子を示す($a=0.004$ はノートと同じ値だが,級数の打ち切り項数は,ノートでは $n=1$〜$100$ までであるところ,本章ではより滑らかな図にするため $n=1$〜$200$ まで計算した).
例題43.6 初期条件が単振動1つで済む場合(熱伝導方程式)
長さ $L$,両端固定の棒が,初期条件 $U(x,0)=T_0\sin\dfrac{\pi x}L$ から出発するとき,熱伝導方程式の解を求めよ.例題43.4(波動方程式,同じ初期条件の形)の結果と比較せよ.
解答 例題43.4と全く同じ直交性の議論により,$b_1=T_0$,$n\ne1$ では $b_n=0$ となる.したがって解は
$$ U(x,t)=T_0\sin\frac{\pi x}Le^{-\frac{\pi^2}{L^2}at} $$である.例題43.4の波動方程式の解 $u(x,t)=h\sin(\pi x/L)\cos(\pi vt/L)$ と,形($\sin(\pi x/L)$ という「山」の形)は全く同じ初期条件から出発しているにもかかわらず,時間発展の仕方は対照的である:波動方程式の解は $\cos(\pi vt/L)$ に従って永久に振動し続けるが,熱伝導方程式の解は $e^{-\pi^2at/L^2}$ に従って時間とともに指数関数的に減衰し,$t\to\infty$ で $U\to0$ に近づく.同じ「空間的な形」を持つ初期条件でも,方程式の種類(時間微分が2階か1階か)によって運命が全く異なるのである——これが公式43.5で述べた質的な違いの,最も単純な具体例である.(sympyで検算:$b_1=T_0$,$n=2,\ldots,5$ で $b_n=0$ となることを直接積分して確認済み.)
43.7 2次元熱伝導方程式の変数分離解法
43.5節で波動方程式を2次元に拡張したのと同じように,熱伝導方程式も2次元(平面上に広がる薄い板の温度分布)に拡張できる.変数分離の手順は「$x,y$ の2つの空間変数についてそれぞれ固有値問題を解き,$t$ については1階の減衰方程式を解く」という,これまでの組み合わせそのものである.
例43.4 2次元熱伝導方程式:帯状の初期高温領域(ノートの例題)
関数 $U(x,y,t)$ について,熱伝導方程式
$$ \pdiff ut=a\left(\pdiff{^2u}{x^2}+\pdiff{^2u}{y^2}\right)\qquad(0\lt x\lt L,\ 0\lt y\lt L,\ 0\lt t) $$を,初期条件
$$ U(x,y,0)=\begin{cases}T_0&\left(\dfrac L2\le x\le\dfrac34L,\ 0\lt y\lt L\right)\\0&\left(0\le x\lt\dfrac L2\ \text{または}\ \dfrac34L\lt x\le L,\ 0\le y\le L\right)\end{cases} $$境界条件 $U(0,y,t)=U(L,y,t)=U(x,0,t)=U(x,L,t)=0$(正方形の板の縁が $0\,{}^\circ$ に保たれている)のもとで解け.
解答(1)変数分離と固有値問題 $U(x,y,t)=X(x)Y(y)T(t)$ とおく.方程式に代入すると $XY\dot T=a(X''YT+XY''T)$,両辺を $aXYT$ で割ると
$$ \frac{\dot T}{aT}=\frac{X''}X+\frac{Y''}Y=-\omega^2\qquad(\omega\gt0) $$(右辺が $x,y$ の関数の和として恒等的に定数になるという理屈は43.4.1項(43.5節でも用いた)と同じ.左辺が正になると $T$ が指数関数的に発散してしまい物理的でないので,定数を $-\omega^2\ (\lt0)$ の形で置くのが自然である——実際,境界条件のもとで非自明解が存在するのはこの符号のときだけであることが,以下の計算でも確認できる).ここでさらに $\dfrac{X''}X=-\omega_1^2$,$\dfrac{Y''}Y=-\omega_2^2$ とおくと,43.5節と全く同じ固有値問題になる.境界条件 $X(0)=X(L)=0$ から
$$ X(x)=B_1\sin\frac{m\pi}Lx\qquad(\omega_1L=m\pi,\ m=1,2,\ldots) $$境界条件 $Y(0)=Y(L)=0$ から
$$ Y(y)=B_2\sin\frac{n\pi}Ly\qquad(\omega_2L=n\pi,\ n=1,2,\ldots) $$が得られる.最後に $T$ の方程式 $\dot T=-a\omega^2T$(43.6節と同じ1階減衰方程式)を解くと
$$ T(t)=A_3e^{-a\omega^2t} $$ここで $\omega^2=\omega_1^2+\omega_2^2=\dfrac{m^2+n^2}{L^2}\pi^2$ であるから,
$$ T(t)=A_3\,e^{-\frac{m^2+n^2}{L^2}a\pi^2t} $$解答(2)重ね合わせとフーリエ係数 以上をまとめると,各 $m,n$ に対する解は
$$ U(x,y,t)=B_1B_2A_3\,e^{-\frac{m^2+n^2}{L^2}a\pi^2t}\sin\frac{m\pi}Lx\sin\frac{n\pi}Ly $$である($B_1B_2A_3=b_{mn}$ とおく).重ね合わせの原理(定理43.1)により,一般解は
$$ U(x,y,t)=\sum_{m=1}^\infty\sum_{n=1}^\infty b_{mn}\,e^{-\frac{m^2+n^2}{L^2}a\pi^2t}\sin\frac{m\pi}Lx\sin\frac{n\pi}Ly $$$t=0$ を代入し,2重フーリエ正弦級数の係数公式を使うと,
$$ b_{mn}=\frac4{L^2}\int_0^L\int_0^LU(x,y,0)\sin\frac{m\pi}Lx\sin\frac{n\pi}Ly\,dx\,dy $$初期条件は $x\in[L/2,3L/4]$ でだけ $T_0$,$y$ 方向には全域で一様なので,2重積分は $x$ 方向の積分と $y$ 方向の積分に分解できる.
$$ b_{mn}=\frac4{L^2}T_0\int_{L/2}^{3L/4}\sin\frac{m\pi}Lx\,dx\int_0^L\sin\frac{n\pi}Ly\,dy $$2つの積分をそれぞれ計算すると,
$$ \int_{L/2}^{3L/4}\sin\frac{m\pi}Lx\,dx=\left[-\frac L{m\pi}\cos\frac{m\pi}Lx\right]_{L/2}^{3L/4}=\frac L{m\pi}\left(\cos\frac{m\pi}2-\cos\frac34m\pi\right) $$ $$ \int_0^L\sin\frac{n\pi}Ly\,dy=\left[-\frac L{n\pi}\cos\frac{n\pi}Ly\right]_0^L=\frac L{n\pi}(1-\cos n\pi)=\frac L{n\pi}\bigl(1-(-1)^n\bigr) $$であるから,
$$ b_{mn}=\frac{4T_0}{mn\pi^2}\left(\cos\frac{m\pi}2-\cos\frac34m\pi\right)\bigl(1-(-1)^n\bigr) $$が得られ,求める解は
\begin{equation} U(x,y,t)=\sum_{m=1}^\infty\sum_{n=1}^\infty\frac{4T_0}{mn\pi^2}\left(\cos\frac{m\pi}2-\cos\frac34m\pi\right)\bigl(1-(-1)^n\bigr)\cdot e^{-\frac{m^2+n^2}{L^2}\pi^2at}\sin\frac{m\pi}Lx\sin\frac{n\pi}Ly \label{eq:43-heat2d-sol} \end{equation}である(例43.4の解答終わり).係数 $(1-(-1)^n)$ は $n$ が奇数のとき $2$,偶数のとき $0$ になることに注目すると,$y$ 方向には奇数番目の固有振動だけが現れることが分かる——これは,初期条件が $y$ 方向には一定(帯全体で同じ高さ)であり,$y=L/2$ を軸とした左右対称な分布になっているためである.
例題43.7 初期条件が単振動1つで済む場合(2次元熱伝導方程式)
正方形 $0\lt x\lt L,\ 0\lt y\lt L$ の板が,初期条件 $U(x,y,0)=T_0\sin\dfrac{\pi x}L\sin\dfrac{2\pi y}L$ から出発するとき,熱伝導方程式の解を求めよ.
解答 初期条件は $m=1,n=2$ の固有関数そのものなので,2重フーリエ正弦級数の直交性から $b_{12}=T_0$,それ以外の $b_{mn}=0$ である.$m^2+n^2=1^2+2^2=5$ に注意すると,求める解は
$$ U(x,y,t)=T_0\sin\frac{\pi x}L\sin\frac{2\pi y}Le^{-\frac{5\pi^2}{L^2}at} $$である.$m,n$ の値が大きいほど($m^2+n^2$ が大きいほど)減衰が速いことに注意しよう——これは「細かい模様(波長の短い温度分布)ほど早くならされる」という直感と一致する.たとえば $(m,n)=(1,1)$ の減衰率は $2a\pi^2/L^2$,$(m,n)=(1,2)$ は $5a\pi^2/L^2$ であり,後者の方が $2.5$ 倍速く減衰する.(sympyで検算:$m=1,n=2$ で $b_{12}=T_0$,それ以外の $(m,n)\,(\le3)$ ではすべて $0$ となることを直接積分して確認済み.)
43.8 ラプラス方程式の例題
本章の最後に,時間 $t$ を全く含まない,もう1つの重要な2階線形偏微分方程式——ラプラス方程式(Laplace's equation)$\pdiff{^2u}{x^2}+\pdiff{^2u}{y^2}=0$ を扱う.これは,熱伝導方程式 $\pdiff ut=a\left(\pdiff{^2u}{x^2}+\pdiff{^2u}{y^2}\right)$ で,十分に時間が経って定常状態($\pdiff ut=0$)に達した場合の方程式でもある(43.3節の1次元の定常解の議論を2次元に一般化したものと考えてもよい).ラプラス方程式は,このあと第44章で座標変換とともにさらに詳しく学び,静電ポテンシャルや流体力学など,物理学のあらゆる場面で登場する.
例43.5 2次元ラプラス方程式の境界値問題(ノートの例題)
関数 $u(x,y)$ について,ラプラス方程式
$$ \pdiff{^2u}{x^2}+\pdiff{^2u}{y^2}=0\qquad(0\lt x\lt L,\ 0\lt y\lt L) $$を,境界条件
$$ u(0,y)=\begin{cases}T_0&\left(0\le y\le\dfrac L2\right)\\0&\left(\dfrac L2\lt y\le L\right)\end{cases}\, ,\qquad u(x,0)=u(x,L)=u(L,y)=0 $$のもとで解け.
解答(1)変数分離 $u(x,y)=X(x)Y(y)$ とおく.ラプラス方程式に代入すると $X''Y+XY''=0$,両辺を $XY\,(\ne0)$ で割ると
$$ \frac{X''}X+\frac{Y''}Y=0\qquad\therefore\ -\frac{X''}X=\frac{Y''}Y=C\qquad(C:\text{定数}) $$境界条件のうち,$y$ について同次(右辺が $0$)なのは $u(x,0)=u(x,L)=0$,すなわち $Y(0)=Y(L)=0$ である($x=0$ での境界条件は同次でないので,これは後で使う).そこでまず,$Y(0)=Y(L)=0$ のもとで $Y''=CY$ の非自明解を探す.
[i] $C\gt0$ のとき,一般解は $Y(y)=A_1e^{\sqrt Cy}+B_1e^{-\sqrt Cy}$ であり,43.4.2項(i)と全く同じ議論により,境界条件 $Y(0)=Y(L)=0$ を満たす非自明な $A_1,B_1$ は存在しない.
[ii] $C=0$ のとき,一般解は $Y(y)=A_2y+B_2$ であり,43.4.2項(ii)と同じ議論により,やはり非自明解は存在しない.
[iii] $C=-\omega^2$($\omega\gt0$)とおくと,43.4.2項(iii)と同じ議論により,$Y(0)=Y(L)=0$ を満たす非自明解が存在するのは $\omega L=n\pi$($n=1,2,\ldots$)のときであり,
$$ Y(y)=B_1\sin\frac{n\pi}Ly $$である.
解答(2)$X(x)$ の決定——ノートの続き $C=-\omega^2$($\omega=n\pi/L$)が定まったので,$X$ の方程式に代入する.
$$ -\frac{X''}X=-\omega^2\qquad\therefore\ X''=\omega^2X $$一般解は
$$ X(x)=A_3e^{\omega x}+B_3e^{-\omega x} $$ここで使う境界条件は,$x=L$ での同次条件 $u(L,y)=0$,すなわち $X(L)=0$ である($x=0$ での境界条件は同次でないので,$X(0)$ には制約を課さない——実際,$u(0,y)$ が $0$ でない値を取れるのは,$X(0)\ne0$ だからこそである).
$$ X(L)=A_3e^{\omega L}+B_3e^{-\omega L}=0\qquad\therefore\ B_3=-A_3e^{2\omega L} $$これを $X(x)$ の式に代入すると,
$$ X(x)=A_3e^{\omega x}-A_3e^{2\omega L}e^{-\omega x}=A_3e^{\omega L}\bigl(e^{\omega(x-L)}-e^{-\omega(x-L)}\bigr) $$(ノートはここまでの計算の直前,「$X(x)=A_3(\ldots)$」の途中で中断している.以下は本章で補って完成させた続きである.)ここで,双曲線関数(第31章のラプラス変換で登場した $\sinh,\cosh$,あるいは第2章で扱った双曲線関数)の定義 $\sinh\theta=\dfrac{e^\theta-e^{-\theta}}2$ を使うと,かっこの中は $e^{\omega(x-L)}-e^{-\omega(x-L)}=2\sinh\bigl(\omega(x-L)\bigr)=-2\sinh\bigl(\omega(L-x)\bigr)$($\sinh$ は奇関数だから符号を反転できる)と書ける.したがって,定数を整理し直すと,
\begin{equation} X(x)=B_4\sinh\bigl(\omega(L-x)\bigr) \label{eq:43-laplace-X} \end{equation}という形にまとめられる($B_4$ は $A_3$ を含む新しい任意定数).この形は,$X(L)=B_4\sinh0=0$ を自動的に満たしており,境界条件 $u(L,y)=0$ が最初から組み込まれていることが分かる(双曲線正弦 $\sinh$ を使うと,指数関数のままより見通しがよくなる).
解答(3)重ね合わせと非同次境界条件の適用 以上から,各 $n=1,2,\ldots$ に対する解は
$$ u_n(x,y)=B_1B_4\sinh\bigl(\omega(L-x)\bigr)\sin\frac{n\pi}Ly\qquad\left(\omega=\frac{n\pi}L\right) $$である($B_1B_4=b_n$ とおく).ラプラス方程式もまた線形なので,重ね合わせの原理により,一般解は
$$ u(x,y)=\sum_{n=1}^\infty b_n\sinh\!\left(\frac{n\pi}L(L-x)\right)\sin\frac{n\pi}Ly $$となる.残る条件は,まだ使っていない $x=0$ での境界条件 $u(0,y)=(\text{段差関数})$ である.$x=0$ を代入すると,$\sinh(n\pi)$ は $x$ に依存しない定数として係数に吸収されて,
$$ u(0,y)=\sum_{n=1}^\infty b_n\sinh(n\pi)\sin\frac{n\pi}Ly=u(0,y) $$となる.これはフーリエ正弦級数の形をしているから,係数 $b_n\sinh(n\pi)$ が,通常のフーリエ正弦係数の公式で決まる.
$$ b_n\sinh(n\pi)=\frac2L\int_0^Lu(0,y)\sin\frac{n\pi}Ly\,dy=\frac2L\int_0^{L/2}T_0\sin\frac{n\pi}Ly\,dy $$右辺の積分を計算すると,
$$ \frac2L\int_0^{L/2}T_0\sin\frac{n\pi}Ly\,dy=\frac2L\left[-\frac L{n\pi}T_0\cos\frac{n\pi}Ly\right]_0^{L/2}=\frac{2T_0}{n\pi}\left(1-\cos\frac{n\pi}2\right) $$したがって $b_n=\dfrac{2T_0}{n\pi\sinh(n\pi)}\left(1-\cos\dfrac{n\pi}2\right)$ であり,求める解は
\begin{equation} u(x,y)=\sum_{n=1}^\infty\frac{2T_0}{n\pi\sinh(n\pi)}\left(1-\cos\frac{n\pi}2\right)\sinh\!\left(\frac{n\pi}L(L-x)\right)\sin\frac{n\pi}Ly \label{eq:43-laplace-sol} \end{equation}である(例43.5の解答終わり,ノートで中断していた計算を完成させた).
なぜ?:$\sinh(n\pi)$ が分母にあるのは自然なこと
$\sinh(n\pi)$ は $n$ が大きくなると指数関数的に急速に大きくなる($\sinh\theta\approx e^\theta/2$,$\theta$ が大きいとき).したがって係数 $b_n$ は $n$ が大きいほど急激に小さくなり,級数はきわめて速く収束する.これは,ラプラス方程式の解が境界のなめらかでない部分($y=L/2$ における段差)の影響を,領域の内部($x$ が $0$ から離れるほど)では急速に「忘れて」なめらかになっていく,という性質を表している——熱伝導方程式の定常状態が,境界の凹凸を内部でなめらかに均してしまうのと同じ現象である.
例題43.8 境界条件を単純化した場合
例43.5と同じ領域・同じ方程式で,境界条件を $u(0,y)=T_0$($0\le y\le L$ の全域で一定),$u(x,0)=u(x,L)=u(L,y)=0$ に単純化したとき,解を求めよ.
解答 解答(1)(2)の議論は境界条件によらないので,一般解の形 $u(x,y)=\displaystyle\sum_nb_n\sinh\!\left(\frac{n\pi}L(L-x)\right)\sin\frac{n\pi}Ly$ はそのまま使える.$x=0$ での境界条件だけが変わるので,フーリエ係数を計算し直す.
$$ b_n\sinh(n\pi)=\frac2L\int_0^LT_0\sin\frac{n\pi}Ly\,dy=\frac2L\left[-\frac L{n\pi}T_0\cos\frac{n\pi}Ly\right]_0^L=\frac{2T_0}{n\pi}\bigl(1-\cos n\pi\bigr)=\frac{2T_0}{n\pi}\bigl(1-(-1)^n\bigr) $$$1-(-1)^n$ は $n$ が奇数のとき $2$,偶数のとき $0$ になるから,奇数番目の項だけが残る.したがって $b_n=\dfrac{4T_0}{n\pi\sinh(n\pi)}$($n$ が奇数のときのみ)であり,求める解は
$$ u(x,y)=\sum_{n\,\text{奇数}}\frac{4T_0}{n\pi\sinh(n\pi)}\sinh\!\left(\frac{n\pi}L(L-x)\right)\sin\frac{n\pi}Ly $$である.境界条件が $y$ 方向に一定($y\to L-y$ の反転に対して対称)になったことで,フーリエ係数が奇数次の項だけに簡単化された——これは例43.4(2次元熱伝導方程式)で見た「$y$ 方向の対称性が奇数次だけを残す」のと同じ現象である.(sympyで検算:$(2/L)\int_0^LT_0\sin(n\pi y/L)dy=2T_0(1-(-1)^n)/(n\pi)$ となることを確認済み.また,級数の各項 $\sinh(n\pi(L-x)/L)\sin(n\pi y/L)$ が $\nabla^2u=0$ を満たすことも直接微分して確認済み.)
43.9 まとめと演習
43.9.1 まとめ
- 真空中のMaxwell方程式(ベクトル恒等式 $\nabla\times(\nabla\times\bm E)=-\nabla^2\bm E+\nabla(\nabla\cdot\bm E)$ を使う)から,電場の波動方程式 $\nabla^2\bm E=\dfrac1{c^2}\pdiff{^2\bm E}{t^2}$(公式43.1)が導かれる.
- 波動方程式に単振動解を仮定し,ド・ブロイの関係式を組み合わせると,時間に依存しないシュレーディンガー方程式(公式43.2)が現れる——第45章への伏線.
- 弦の微小要素の運動方程式 $ma=F$ から,弦の波動方程式 $\pdiff{^2u}{x^2}=\dfrac1{v^2}\pdiff{^2u}{t^2}$($v^2=T/\rho$,公式43.3)が導かれる.
- 棒の熱収支とフーリエの熱伝導法則から,熱伝導方程式 $\pdiff ut=a\pdiff{^2u}{x^2}$($a=k/(\rho C)$,公式43.4)が導かれる.
- 変数分離法(定義43.1)$u=X(x)T(t)$(あるいは $X(x)Y(y)T(t)$)を偏微分方程式に代入すると,分離定数を通じていくつかの常微分方程式に分解される.境界条件のもとで非自明解が存在するのは,分離定数がとびとびの値(固有値)を取るときだけであり,そのときの解を固有関数と呼ぶ.
- 重ね合わせの原理(定理43.1):線形偏微分方程式では,固有関数の無限1次結合もまた解になる.初期条件・境界条件から,フーリエ級数展開によって係数を決定する.
- 波動方程式では $T(t)=\cos(\cdots)$(永久に振動),熱伝導方程式では $T(t)=e^{-(\cdots)t}$(指数的に減衰)——同じ空間方向の固有関数を持ちながら,時間発展の性質は対照的である(公式43.5).
- 1次元・2次元の波動方程式,1次元・2次元の熱伝導方程式を,具体的な初期条件(三角形分布,山型曲面,矩形パルス,帯状分布)のもとで実際に解き切った.
- 2次元ラプラス方程式の境界値問題を,双曲線関数 $\sinh$ を使って解いた(ノートで中断していた計算を本章で完成させた).
43.9.2 演習問題
演習43.1 波動方程式の速さの次元
密度 $\rho$(線密度),張力 $T$ の弦の波動方程式の速さ $v$ の次元を,$v=\sqrt{T/\rho}$ の右辺の次元を調べることで確認せよ.
ヒント:$T$ の単位は $\mathrm N=\mathrm{kg\cdot m/s^2}$,$\rho$ の単位は $\mathrm{kg/m}$ である.
演習43.2 熱伝導方程式の温度伝導率の次元
熱伝導率 $k$(単位 $\mathrm{W/(m\cdot K)}=\mathrm{kg\cdot m/(s^3\cdot K)}$),密度 $\rho$(単位 $\mathrm{kg/m^3}$),比熱 $C$(単位 $\mathrm{J/(kg\cdot K)}=\mathrm{m^2/(s^2\cdot K)}$)から,温度伝導率 $a=k/(\rho C)$ の次元を求めよ.
ヒント:それぞれの単位を分数の形で代入し,約分する.
演習43.3 波動方程式:初期条件が2倍振動数の正弦波の場合
長さ $L$,両端固定の弦が,初期条件 $u(x,0)=h\sin\dfrac{2\pi x}L$,$\pdiff ut(x,0)=0$ から出発するとき,波動方程式の解を求めよ.
ヒント:例題43.4と同じ考え方が使える.初期条件はどの固有関数($m$ の値)に対応するか.
演習43.4 熱伝導方程式:初期温度が一様な場合
長さ $L$,両端固定($U(0,t)=U(L,t)=0$)の棒が,初期条件 $U(x,0)=T_0$($0\le x\le L$ の全域で一定)から出発するとき,熱伝導方程式の解を求めよ.
ヒント:一般解の形は43.6節と同じ.フーリエ正弦級数の係数 $b_n=\dfrac2L\displaystyle\int_0^LT_0\sin\dfrac{n\pi}Lx\,dx$ を計算せよ.一様な初期条件でも,境界で $U=0$ を強制しているので,境界付近では初期条件との食い違いが生じることに注意(ギブス現象の原因).
演習43.5 2次元波動方程式:異なる固有振動の場合
単位正方形 $0\lt x\lt1,\ 0\lt y\lt1$ の膜が,初期条件 $u(x,y,0)=\sin\pi x\sin2\pi y$,$\pdiff ut(x,y,0)=0$ から出発するとき,波動方程式の解を求めよ.
ヒント:例題43.5と同じ考え方だが,今回は $m=1,n=2$ の固有関数である.
演習43.6 ラプラス方程式:境界条件を2次関数にした場合
例43.5と同じ領域・同じ方程式で,境界条件を $u(0,y)=y(L-y)$($0\le y\le L$),$u(x,0)=u(x,L)=u(L,y)=0$ にしたとき,解を求めよ.
ヒント:一般解の形はそのまま使える.フーリエ係数 $b_n\sinh(n\pi)=\dfrac2L\displaystyle\int_0^Ly(L-y)\sin\dfrac{n\pi}Ly\,dy$ を,部分積分を2回繰り返して計算せよ.
43.9.3 参考文献
- 望月泰英『数学ノート 偏微分方程式』(手書き講義ノート).本章の底本.
- 寺沢寛一『自然科学者のための数学概論』岩波書店.
- 溝畑茂『偏微分方程式論』岩波書店.
- 砂川重信『理論電磁気学』紀伊國屋書店.