1
波動方程式（1 次元空間）
u = u(x, t)に対し utt= a2uxx, a > 0 (1) の形で表される方程式を波動方程式と言う。 波動方程式は、弦や膜などの振動現象を表す。 考える領域は D :={(x, t) : 0 ≤ x ≤ l, t ≥ 0} とする。境界条件は u (0, t) = 0, u (l, t) = 0, t≥ 0 (2) とし、初期条件は u(x, 0) = f (x), ∂u(x, 0) ∂t = 0, 0≤ x ≤ l (3) とする。2
変数分離法による解法
変数分離法では、はじめに u(x, t) = X(x)T (t)の形の u で(1)を満たす ものを求める。代入して整理すると 1 a2 T′′(t) T (t) = X′′(x) X(x) が得られる。これがある定数α に等しいとす ると、二つの常微分方程式 T′′(t) = αa2T (t) (4) X′′(x) = αX(x) (5) が得られる。境界条件 (2)より X(0) = 0, X(l) = 0 (6) が成立する。また初期条件(3) より T (0) = f (x), T′(0) = 0 (7) が成立する。条件(6), (7) の下で、常微分方 程式 (4)と (5)を解くことが、第一の目標と なる。 (5)を満たす X(x)は X(x) = Ae√αx+ Be−√αx (α > 0) Ax + B (α = 0) A cos√−αx + B sin√−αx (α < 0) となる。但し A, B は定数である。境界条件 (6)を満たすもので、自明な解（X(x)≡ 0） でないものは X(x) = Bnsin nπx l ( α =−n 2π2 l2 ) のみである。但しn = 1, 2,· · · , nで、Bn は 定数である。一方、T (t) については T (t) = Cncos nπat l (Cn:定数) が分かる。これらの積を取って u(x, t) = BnCnsin nπx l cos nπat l が、(1) を満たす（一つの）u であることが分 かる。3
重ね合わせの原理
線型の微分方程式(1) を満たすu を足し上げ たものも(1) を満たす。よって u(x, t) = ∞ ∑ n=1 Ansin nπx l cos nπat l も(1)を満たす。但し An:= BnCn とした。4
フーリエ展開の原理
初期条件(7) より f (x) = ∞ ∑ n=1 Ansin nπx l が得られる。フーリエ展開の原理から An= 2 l ∫ l 0 f (x) sinnπx l dx が得られる。以上をまとめて u(x, t) = 2 l ∞ ∑ n=1 ∫ l 0f (y) sinnπy l dy sin nπx l cos nπat l が、条件(2)-(3)の下での波動方程式 (1)の解 となる。
「偏微分方程式」プリント（第
4
回・
2015
年
5
月
11
日）
担当教員：國谷紀良 研究室：工学部本館3W-406 連絡先：tkuniya@port.kobe-u.ac.jp1
波動方程式（2 次元空間）
xy 平面の近くで振動する膜を考える。時刻 t における膜上の点P の座標をz(x, y, t) とし たとき、膜の振動は次の波動方程式で表さ れる。 ∂2z ∂t2 = c 2 ( ∂2z ∂x2 + ∂2z ∂y2 ) , c > 0. (1) 考える領域は D :原点を中心とする半径 aの円 とする。この場合、円柱座標x(r, θ) = r cos θ, y(r, θ) = r sin θ z(x, y, t) = z(r, θ, t) (0≤ r ≤ a, 0 ≤ θ ≤ 2π) を利用すると便利である。このとき(1) は次 のように変換される。 ∂2z ∂t2 = c 2 ( ∂2z ∂r2 + 1 r ∂z ∂r + 1 r2 ∂2z ∂θ2 ) (2) 演習問題1 (1)は円柱座標の下で (2)に変換されるこ とを示せ。 (2)の初期条件は z(r, θ, 0) = φ(r), ∂z(r, θ, 0) ∂t = 0 (3) とし、境界条件は z (a, θ, t) = 0 (4) とする。条件(3)-(4) の下で(2)を解くことを 考える。
2
変数分離法による解法
初期条件がθ に依存しないため z(r, θ, t) = z(r, t) と見なすことが出来る。よって変数分離法に 基づき z(r, t) = R(r)T (t) (5) という形の解を探す。これを(2) に代入して 整理すると T′′(t)− c2kT (t) = 0 R′′(r) + 1 rR ′(r)− kR(r) = 0 (6) （k：定数）という形の二つの常微分方程式を 導出することが出来る。 演習問題2 (5)を(2)に代入して整理することで、(6) を導け。 微分方程式(6)の一般解を求めて、初期条件 (3)と境界条件(4) を利用することで T (t) = c1cos cαn a t と R(r) = J0 (α n a r ) が得られる。但しJ0(·)は位数 0 の第一種 ベッセル関数で、αn はJ0(αn) = 0 を満たす 無数の正根 α1< α2 <· · · < αn <· · · とする。よって定数An に対し z(r, t) = ∞ ∑ n=1 AnJ0 (α n a r ) coscαn a t が求める解となるが、An を決定する必要があ る。初期条件 (3)より φ(r) = ∞ ∑ n=1 AnJ0 (αn a r ) となる。ベッセル関数の公式を用いると An= 2 {J1(αn)}2 ∫ 1 0 sφ (αs) J0(αs) ds が得られる。1
マクスウェルの電磁場方程式
真空中の電場ベクトルE、磁場ベクトル Hに 対するマクスウェルの電磁場方程式は ϵ0 ∂E ∂t = rot H− J0, ϵ0 div E = ρ0 µ0 ∂H ∂t =−rot E − Jm, µ0 div H = ρm (1) となる。ただし ρ0：電荷密度 ρm：磁場密度 J0：電流密度 Jm：磁流密度 であり、ϵ0, µ0 は定数で 1 √ϵ 0µ0 = c = 2.99792× 108 m/s = 光速度 を満たすものとする。経験的に ρm= 0, Jm = 0 (2) と考えられている。物質があるときは、電気 的分極 Pと磁気的分極Mの影響で ρ0 → ρ0− div P J0 → J0+ ∂P ∂t + rotM H→ H + M (3) という置換を施す必要がある。さらに電束密 度D と磁束密度B を D = ϵ0E + P, B = µ0(H + M) (4) で定義すると、(1) は ∂D ∂t = rot H− J0, div D = ρ0 ∂B ∂t =−rot E, µ0 div B = 0 (5) と変換される。マクスウェルの方程式として はこちらの形の方がより知られている。 演習問題1 (2)-(4) の下で、(1)を(5)に変換せよ。2
波動方程式の導出
電磁場は、スカラーポテンシャル V とベクト ルポテンシャル AによってB = rot A, E =−grad V −∂A
∂t (6) のように表すことが出来る。また条件 div A + ϵ0µ0 ∂V ∂t = 0 (7) を課すと、(1)の上段の式は次の波動方程式に 書き換えられる： ∂2V ∂t2 = 1 ϵ0µ0 ∆V + ρ0 ϵ2 0µ0 ∂2A ∂t2 = 1 ϵ0µ0 ∆A + 1 ϵ0 J0 (8) またダランベールの演算子 □ := ∆ − ϵ0µ0 ∂2 ∂t2 を用いると、(8) は □V = −ρ0 ϵ0 , □A = −µ0J0 のように書ける。 演習問題2 (6)-(7)の下で、(1) を(8)に変換せよ。
補足（記号の復習など）
• ベクトル E = (E1, E2, E3) に対し rot E = ( ∂E3 ∂y − ∂E2 ∂z , ∂E1 ∂z − ∂E3 ∂x , ∂E2 ∂x − ∂E1 ∂y ) div E = ∂E1 ∂x + ∂E2 ∂y + ∂E3 ∂z • スカラー関数 V = V (x, y, z)に対し grad V = ( ∂V ∂x, ∂V ∂y, ∂V ∂z ) • スカラー V に対し div (grad V ) = ∆V • ベクトル A = (A1, A2, A3)に対し「偏微分方程式」プリント（第
6
回・
2015
年
5
月
25
日）
担当教員：國谷紀良 研究室：工学部本館3W-406 連絡先：tkuniya@port.kobe-u.ac.jp1
波動方程式（3 次元空間）
電磁波や音波など、3次元的に広がる波動量 u(x, y, z, t) に対しては、次の3次元波動方程 式を考えることが出来る。 ∂2u ∂t2 = c 2 ( ∂2u ∂x2 + ∂2u ∂y2 + ∂2u ∂z2 ) , c > 0. (1)2
球面波
一点からすべての方向に同じ速さで広がる波 を球面波という。このとき uは、原点からの 距離 r =√x2+ y2+ z2 と 時間t の関数 u = u (r, t) となる。計算すると ∂2u ∂x2 = 1 r ∂u ∂r − x2 r3 ∂u ∂r + x2 r2 ∂2u ∂r2 (2) ∂2u ∂y2 = 1 r ∂u ∂r − y2 r3 ∂u ∂r + y2 r2 ∂2u ∂r2 (3) ∂2u ∂z2 = 1 r ∂u ∂r − z2 r3 ∂u ∂r + z2 r2 ∂2u ∂r2 (4) となる。r2= x2+ y2+ z2 より ∂2u ∂x2 + ∂2u ∂y2 + ∂2u ∂z2 = 1 r ∂2(ru) ∂r2 (5) が得られる。 演習問題1 (2)-(4) より (5)を求め、(1)は ∂2(ru) ∂t2 = c 2∂2(ru) ∂r2 (6) に書き換えられることを示せ。 (6)は1次元波動方程式と見なせる。その一般 解は、任意関数 g1, g2 に対し u (r, t) = 1 rg1(r− ct) + 1 rg2(r + ct) (7) となる。この右辺第一項は速さ cで波元から 広がる発散波を、第二項は波元に集まる収束 波を表す。 演習問題2 (7)で与えられるu(r, t)は1次元波動方程 式(6)を満たすことを示せ。3
平面波
一方向にある速さで平面として進む波を平面 波という。単位ベクトル (l, m, n)の方向に速 さc で進む(1) の平面波は u (x, y, z, t) = f (lx + my + nz∓ ct) と書ける。例えば、定数 A, k に対する正弦波 u (x, y, z, t) = A sin k (lx + my + nz∓ ct) (8) は(1)を満たす。A を振幅、 k := (kl, km, kn) を波数ベクトル、k (lx + my + nz∓ ct)を波 の位相という。特にr = (x, y, z) を位置ベク トルとすると、(8)は u = A sin (k· r ∓ kct) と書き換えられる。 演習問題3 定数 A, k に対する余弦波 u (x, y, z, t) = A cos k (lx + my + nz∓ ct) = A cos (k· r ∓ kct) は(1)を満たすことを示せ。 オイラーの公式 eix = cos x + i sin xに注意す ると、平面波は一般に u = A exp (i (k· r ∓ kct)) と表示される。2
次元円板におけるディリクレ問題
u = u(x, y)に対し、 uxx+ uyy = 0 (あるいは∆u = 0) (1) をラプラス方程式という。ラプラス方程式は、 時間不変な分布（温度や静電場など）を表す。 u は閉円板 U ={(x, y)| x2+ y2 ≤ 1} 上で連続とし、境界x2+ y2 = 1 において u(x, y) = F (x, y) (2) を満たすものとする。(2) のような条件をディ リクレ境界条件といい、対応する問題(1)-(2) をディリクレ問題という。 極座標変換 x = r cos θ, y = r sin θ によってu = u(r, θ) を考える。このとき uxx+ uyy = urr+ 1 rur+ 1 r2uθθ に注意する（第4回のプリントを参照）と、 問題 (1)-(2) は r2urr+ rur+ uθθ= 0, 0≤ r < 1 u(1, θ) = f (θ), 0≤ θ ≤ 2π u(r, θ) = u(r, θ + 2π) (3) と書き直せる。但し境界 x2+ y2= r2 = 1 上で F (x, y) = F (r, θ) = F (1, θ) = f (θ) とした。このとき、変数分離法に基づき u(r, θ) = R(r)Θ(θ) を(3) に代入して整理することで、次の二つ の常微分方程式が得られる： Θ′′(θ) + αΘ = 0 r2R′′(r) + rR′(r)− αR(r) = 0 α :定数 (3)の第三式より Θ(θ) = Θ(θ + 2π) が成立す ることに注意すると、Θ(θ)の一般解は Θ(θ) = Ancos nθ + Bnsin nθ, n = 1, 2,· · · となる。ただし α = n2 で、R(r)の一般解は R(r) = Cnrn となる。よって an= AnCn, bn= BnCn とす ればu(r, θ) = rn(ancos nθ + bnsin nθ) , n = 1, 2,· · · が得られる。特に重ね合わせの原理によって u(r, θ) = a0 2 + ∞ ∑ n=1 rn(ancos nθ + bnsin nθ) も方程式の解になる。境界条件 u(1, θ) = f (θ) より f (θ) = a0 2 + ∞ ∑ n=1 (ancos nθ + bnsin nθ) が得られるが、これはフーリエ級数展開なので an= 1 π ∫ 2π 0 f (φ) cos nφdφ bn= 1 π ∫ 2π 0 f (φ) sin nφdφ が得られる。以上をまとめると、形式解 u(r, θ) = 1 2π ∫ 2π 0 f (φ)dφ +1 π ∞ ∑ n=1 rn ( cos nθ ∫ 2π 0 f (φ) cos nφdφ + sin nθ ∫ 2π 0 f (φ) sin nφdφ ) (4) を得ることが出来る。 演習問題 f (θ) = sin2θ の時の、ラプラス方程式の ディリクレ問題の形式解(4) を求めよ。
「偏微分方程式」演習問題解答（第
7
回・
2015
年
6
月
1
日）
担当教員：國谷紀良 研究室：工学部本館3W-406 連絡先：tkuniya@port.kobe-u.ac.jp演習問題
f (θ) = sin2θ の時の、ラプラス方程式のディ リクレ問題の形式解(3) を求めよ。解答
(3)の右辺の積分を先に計算する。 ∫ 2π 0 f (φ)dφ = ∫ 2π 0 sin2φdφ = ∫ 2π 0 1− cos 2φ 2 dφ = [ φ 2 − sin 2φ 4 ]2π 0 = π. ∫ 2π 0 f (φ) cos nφdφ = ∫ 2π 0 1− cos 2φ 2 cos nφdφ = 1 2 ∫ 2π 0 {cos nφ−cos(n + 2)φ + cos(n− 2)φ 2 } dφ = −1 4 ∫ 2π 0 cos(n− 2)φdφ = { −π 2, n = 2, 0, n̸= 2. ∫ 2π 0 f (φ) sin nφdφ = ∫ 2π 0 1− cos 2φ 2 sin nφdφ = 0. よって、これらを(3) の右辺に代入すると u(r, θ) = 1 2ππ + 1 πr 2cos 2θ(−π 2 ) = 1 2 ( 1− r2cos 2θ).
1
デルタ関数
イギリスの理論物理学者であるディラック （1902∼1984）は、次の性質を持つ関数を考案 した。 ∫ ∞ −∞δ(x)dx = 1, δ(x) = 0 (x̸= 0) これをディラックのデルタ関数という。実連続 関数 f (x)に対し、次が成立する。 ∫ ∞ −∞f (x)δ(x− ξ)dx = f(ξ) 中心 0、分散 s2 の正規分布の密度関数 fs(x) = 1 √ 2πsexp ( −x2 2s2 ) (1) に対し、s→ 0 とするとデルタ関数δ(x) が得 られる（図1）。 -0.50 -0.25 0 0.25 0.5 5 10 15 x s=1/10 s=1/20 s=1/30 s=1/40 図1. 密度関数fs(x)のデルタ関数δ(x)への収束 例題1 (1) で与えられる密度関数は ∫ ∞ −∞fs(x)dx = 1 を満たすことを示せ。2
ポアソン方程式
u = u(x, y)に対し、 uxx+ uyy=−ρ(x, y) (2) をポアソン方程式という。ここでρ(x, y) は湧 き出し量を表す。ρ(x, y) ≡ 0 のときは(2) は ラプラス方程式になる。 ポアソン方程式 (2) の解を求める上で、グ リーン関数を利用することが出来る。3
グリーン関数
デルタ関数を用いて、(2) に対応する次の方 程式を考える。 Gxx+ Gyy=−δ(x)δ(y) (3) この G = G (x, y) をグリーン関数という。も しこのグリーン関数が求められれば、ポアソン 方程式 (2)の解は u(x, y) = ∫∫ G (x− ξ, y − η) ρ (ξ, η) dξdη (4) という形で与えられる。 例題2 (4) で与えられる関数u はポアソン方 程式 (2)を満たすことを示せ。 ポアソン方程式 (2) は、位置ベクトル r を 用いて、次のように表すことも出来る。 ∆u(r) =−ρ(r) (5) このとき(3) と(4)はそれぞれ ∆G(r) =−δ(r) u(r) = ∫ G(r− r′)ρ(r′)dr′ と表すことが出来る。フーリエ変換を利用する と、二次元（r = (x, y)）の場合は G(r) =− 1 2πlog|r| (6) であり、三次元（r = (x, y, z)）の場合は G(r) = 1 4π|r| (7) であることが分かる。 例題3 三次元の場合のグリーン関数(7) は、 r̸= 0 に対して ∆G(r) = 0を満たすことを示 せ（デルタ関数の性質）。 演習問題 二次元の場合のグリーン関数(6)は、r̸= 0 に対して∆G(r) = 0を満たすことを示せ。「偏微分方程式」演習問題解答（第
8
回・
2015
年
6
月
8
日）
担当教員：國谷紀良 研究室：工学部本館3W-406 連絡先：tkuniya@port.kobe-u.ac.jp 演習問題 二次元の場合のグリーン関数(6)は、r̸= 0 に対して∆G(r) = 0を満たすことを示せ。解答
r = (x, y)とすると、(6)は G(r) = G(x, y) =− 1 2πlog √ x2+ y2 と表すことが出来る。よって Gx(x, y) = − 1 2π√x2+ y2 2x 2√x2+ y2 = − x 2π(x2+ y2) Gy(x, y) = − 1 2π√x2+ y2 2y 2√x2+ y2 = y 2π(x2+ y2) となるので、 Gxx(x, y) = − x2+ y2− x · 2x 2π(x2+ y2)2 = x 2− y2 2π(x2+ y2)2 Gyy(x, y) = − x2+ y2− y · 2y 2π(x2+ y2)2 = y 2− x2 2π(x2+ y2)2. よって ∆G(r) = Gxx+ Gyy = 0 (r̸= 0).1
フーリエ正弦級数展開
区間[0, a]上で定義された区分的に連続な関 数f (x)を、奇関数として全実数区間に拡張す るとき、次の級数展開が可能となる。 f (x) = ∞ ∑ n=1 bnsin nπ a x ただし bn= 2 a ∫ a 0 f (x) sinnπ a x dx, n = 1, 2,· · · である。この展開をフーリエ正弦級数展開と いう。 例題1 関数 f (x) = cos x (0≤ x ≤ π)をフー リエ正弦級数に展開せよ。2
長方形領域のラプラス方程式
長方形領域 D ={(x, y) : 0 ≤ x ≤ a, 0 ≤ y ≤ b} (1) を考える。Dの内部では、ラプラス方程式 ∆u = uxx+ uyy = 0 を満たし、境界 ∂D では次の境界条件 u(x, 0) = 0, u(x, b) = 0, 0≤ x ≤ au(0, y) = φ(y), u(a, y) = 0, 0≤ y ≤ b (2) を満たす解u = u(x, y)を求める。変数分離法 を使うと、解の形は u(x, y) = ∞ ∑ n=1 Ansinh nπ b (a− x) sin nπ b y であることが分かる。x = 0とすると、境界条 件(2) より φ(y) = ∞ ∑ n=1 Ansinh nπ b a sin nπ b y であることが分かる。 例題2 フーリエ正弦級数展開を用いて、係数 An を決定せよ。
3
長方形領域のポアソン方程式
長方形領域Dの内部では、ポアソン方程式 ∆u = uxx+ uyy =−ρ(x, y) (3) を満たし、境界 ∂D では次の境界条件 u(x, 0) = 0, u(x, b) = 0, 0≤ x ≤ a u(0, y) = 0, u(a, y) = 0, 0≤ y ≤ b (4) を満たす解 u = u(x, y) を求める。次の形の 関数 u(x, y) = ∞ ∑ m=1 ∞ ∑ n=1 Amnsin mπ a x sin nπ b y を考えると、これは境界条件 (4) を満たすの で、これがポアソン方程式(3)を満たすように 係数 Amn を決定すればよい。 ρ(x, y) をx についてフーリエ正弦級数展開 すると ρ(x, y) = ∞ ∑ m=1 ρm(y) sin mπ a x となる。ただし ρm(y) = 2 a ∫ a 0 ρ(ξ, y) sinmπ a ξ dξ, m = 1, 2,· · · である。さらに ρm(y) をy についてフーリエ 正弦級数展開し、(3) に代入して係数比較をす ると、最終的に Amn= Bmn (mπ a )2 + (nπ b )2 Bmn= 4 ab ∫ a 0 ∫ b 0 ρ(ξ, η) sinmπ a ξ sin nπ b η dηdξ であることが分かる。 演習問題 ρ(x, y) = 1 4xy のとき、境界条件 (4) を満 たすポアソン方程式(3)の解 uを求めよ。「偏微分方程式」演習問題解答（第
9
回・
2015
年
6
月
15
日）
担当教員：國谷紀良 研究室：工学部本館3W-406 連絡先：tkuniya@port.kobe-u.ac.jp 演習問題 ρ(x, y) = 1 4xy のとき、境界条件 (4) を満 たすポアソン方程式(3)の解 uを求めよ。解答
ρ(x, y) = 1 4xy なので Bmn = 1 ab (∫ a 0 ξ sinmπ a ξdξ ) (∫ b 0 η sinnπ b ηdη ) ここで ∫ a 0 ξ sinmπ a ξdξ = [ − a mπξ cos mπ a ξ ]a 0+ a mπ ∫ a 0 cosmπ a ξdξ = −a 2 mπcos mπ + a mπ [ a mπsin mπ a ξ ]a 0 = −a 2 mπ(−1) m = a 2 mπ(−1) m+1 同様に ∫ b 0 η sinnπ b ηdη = b2 nπ(−1) n+1 なので、 Bmn = 1 ab a2 mπ(−1) m+1 b2 nπ(−1) n+1 = ab mnπ2(−1) m+n よって Amn= (−1)m+n (mπ a )2 + (nπ b )2 ab mnπ2 であり、 u(x, y) = ∞ ∑ m=1 ∞ ∑ n=1 (−1)m+n (mπ a )2 + (nπ b )2 ab mnπ2 × sinmπ a x sin nπ b y1
3
次元領域のラプラス方程式
3次元領域のラプラス方程式は
∆u = uxx+ uyy+ uzz = 0 (1) となる。変数分離法に基づき、
u(x, y, z) = X(x)Y (y)Z(z)
を代入して整理すると、三つの常微分方程式 X′′(x) = k12X(x), Y′′(x) = k22Y (y) Z′′(x) = k23Z(z) (2) が得られる。ここでk1, k2, k3 は k21+ k22+ k23 = 0 (3) を満たす複素数である。(2) の一般解は X(x) = c1ek1x+ c′1e−k1x Y (y) = c2ek2y+ c′2e−k2y Z(z) = c3ek3z+ c′3e−k3z となる。ここで c1, c′1, c2, c′2, c3, c′3 は任意定数 である。よって(1)の一般的な解（の一つ）と して u(x, y, z) = ∑ k1,k2,k3 Ck1k2k3e k1x+k2y+k3z が得られる。ここで Ck1k2k3 は任意定数であ り、和 ∑k 1k2k3 は (3) を満たすすべての複素 数k1, k2, k3 に対して取られる。
2
円柱座標に対する変数分離法
円柱座標(ρ, φ, z)は、座標変換 x = ρ cos φ, y = ρ sin φ, z = z (4) で得られる。このとき、ラプラス方程式(1)は uρρ+ 1 ρuρ+ 1 ρ2uφφ+ uzz = 0 (5) と変換される。 例題1 ラプラス方程式 (1)に対する座標変換 (4)の下で (5)を導出せよ。 変数分離法に基づき u (ρ, φ, z) = R(ρ)Φ(φ)Z(z) を(5)に代入して整理すると、微分方程式 R′′(ρ) +R ′(ρ) ρ + ( α2−m 2 ρ2 ) R(ρ) = 0 (6) Φ′′(φ) =−m2Φ(φ), Z′′(z) = α2Z(z) が得られる。ここでα, mは定数である。Φ(φ), Z(z) の一般解は Φ(φ) = A cos mφ + B sin mφ Z(z) = Ceαz+ De−αz (A, B, C, D :定数) となる。R(ρ) の方程式(6) は、α = 1のとき ベッセル微分方程式と呼ばれ、その解は m 次 のベッセル関数と呼ばれる。3
球座標に対する変数分離法
球座標(r, θ, φ) は、座標変換x = r sin θ cos φ, y = r sin θ sin φ z = r cos θ で得られる。このとき、ラプラス方程式(1)は urr+ 2 rur+ 1 r2uθθ+ 1 r2tan θuθ+ 1 r2sin2θuφφ= 0 と変換される。変数分離法に基づき u (r, θ, φ) = R(r)Θ(θ)Φ(φ) を代入して整理すると、一般解
Φ(φ) = A sin mφ + B cos mφ (A, B :定数)
および微分方程式 R′′(r) + 2R ′(r) r − n (n + 1) R(r) r2 = 0 (7) d dµ { (1− µ2)dΘ dµ } + { n(n + 1)− m 2 1− µ2 } Θ = 0 が得られる。ただし µ = cos θである。これら はそれぞれ、オイラーの微分方程式とルジャン ドル陪微分方程式と呼ばれる。 演習問題 オイラーの微分方程式(7)に対し、R(r) = rλ（λ：未定定数）の形を持つ基本解を求 めよ。
「偏微分方程式」演習問題解答（第
10
回・
2015
年
6
月
22
日）
担当教員：國谷紀良 研究室：工学部本館3W-406 連絡先：tkuniya@port.kobe-u.ac.jp 演習問題 オイラーの微分方程式(7)に対し、R(r) = rλ（λ：未定定数）の形を持つ基本解を求 めよ。解答
R(r) = rλ をオイラーの微分方程式に代入す ると λ (λ− 1) rλ−2+ 2λrλ−2− n (n + 1) rλ−2 = 0 より、 λ (λ− 1) + 2λ − n(n + 1) = 0 が得られる。整理すると (λ− n) (λ + n + 1) = 0 より、λ = n,−n − 1が得られる。よって求め る基本解は R(r) = rn, r−n−1 となる。1
有限区間の熱方程式
熱方程式は ut= kuxx, k > 0 (1) で与えられる。ここでu = u(x, t)は時間 t に おける位置 x の温度を表し、k は熱拡散率と 呼ばれる。長さℓの棒における熱の移動を考え る場合、x ∈ [0, ℓ] として(1) は有限区間の熱 方程式となる。 境界条件として u (0, t) = u (ℓ, t) = 0 (2) を課す。これは棒の両端の温度が常に0に保た れていることを意味する。初期条件は u (x, 0) = φ(x) (3) とする。但し φ(x)は連続かつその導関数は区 分的に連続で、φ(0) = φ(ℓ) = 0 を満たすもの とする。 例題1 変数分離法を用いて、条件 (2)-(3) の 下での熱方程式 (1)の解は u (x, t) = ∞ ∑ n=1 2 ℓ ∫ ℓ 0 φ(y) sinnπ ℓ y dy × exp ( −kn2π2 ℓ2 t ) sinnπ ℓ x となることを示せ。 解答の手順（復習） 1) u(x, t) = X(x)T (t) を(1) に代入し、分 離定数を用いて整理することで、X(x)と T (t)それぞれに関する常微分方程式を導 出する。 2) 境界条件 (2) を下に X(x) とT (t) を求 め、その線型和として（任意定数を含む 形で）u(x, t) を表現する。 3) t = 0 のときの関係式から、フーリエ級 数展開を利用して、任意定数を初期条件 φ(x)を含む形で表す。2
無限区間の熱方程式
十分長い棒を想定し、無限区間x∈ (−∞, ∞) での熱方程式(1)を考える。このとき、棒の両 端での温度の影響はないと考え、条件(2)は課 さないことにする。 例題1と同様に変数分離法を利用すると、次 の形の解が得られる。u(x, t) = (a cos λx + b sin λx) e−kλ2t. ここでa, bは任意定数であり、分離定数は−λ2 として計算した。特に a, b はλに依存して連 続的に変化する任意関数a(λ), b(λ)と考えると
uλ(x, t) ={a(λ) cos λx + b(λ) sin λx} e−kλ 2t という特殊解が得られる。このとき、例題1で の線型和に相当して、次の積分も(1)を満たす。 u(x, t) = ∫ ∞ −∞uλ(x, t) dλ = ∫ ∞
−∞{a(λ) cos λx + b(λ) sin λx} e −kλ2t
dλ (4)
初期条件(3) より次が分かる。
φ(x) =
∫ ∞
−∞{a(λ) cos λx + b(λ) sin λx} dλ この式とフーリエ積分公式の係数を比較する と、次が分かる。 a(λ) = 1 2π ∫ ∞ −∞φ(ξ) cos λξdξ (5) b(λ) = 1 2π ∫ ∞ −∞φ(ξ) sin λξdξ (6) これを(4)に代入して整理すると、次のような 無限区間における熱方程式の解が得られる。 u(x, t) = 1 2π ∫ ∞ −∞dλ ∫ ∞ −∞φ(ξ) cos λ (x− ξ) e −kλ2t dξ (7) 演習問題 (5), (6)を (4)に代入して整理することで (7)を求めよ。 ※7月6日（月）の3 限は休講です。
「偏微分方程式」演習問題解答（第
11
回・
2015
年
6
月
29
日）
担当教員：國谷紀良 研究室：工学部本館3W-406 連絡先：tkuniya@port.kobe-u.ac.jp 演習問題 (5), (6)を (4)に代入して整理することで (7)を求めよ。解答
u(x, t) = ∫ ∞ −∞ [{ 1 2π ∫ ∞ −∞φ(ξ) cos λξdξ } cos λx + { 1 2π ∫ ∞ −∞φ(ξ) sin λξdξ } sin λx ] e−kλ2tdλ = 1 2π ∫ ∞ −∞dλ ∫ ∞−∞φ(ξ) (cos λξ cos λx + sin λξ sin λx)
×e−kλ2t dξ = 1 2π ∫ ∞ −∞dλ ∫ ∞ −∞φ(ξ) cos λ (ξ− x) e −kλ2t dξ
1
非斉次な境界条件の熱方程式
有限区間x∈ [0, ℓ]の熱方程式 ut= kuxx, k > 0 (1) を考える。但しu = u(x, t)は点x の時間tで の温度とする。棒の両端x = 0, x = ℓは、それ ぞれ一定温度 u¯0, ¯uℓ の外界にさらされている ものとし、それぞれの熱交換係数を h0, hℓ と すると K ∂u ∂x x=0 = h0( u|x=0− ¯u0) −K ∂u ∂x x=ℓ = hℓ( u|x=ℓ− ¯uℓ) (2) が成り立つ。ここで K は定数（熱伝導率）で ある。 ¯ u0 ̸= 0あるいは u¯ℓ ̸= 0の場合、(2)は非斉 次な境界条件と呼ばれる。変数分離法が適用で きるように、次の関係式 u(x, t) = w(x, t) + α + βx (3) を満たす新たな関数 w = w(x, t) を導入する。 ここでα, βは定数係数である。(3)より、次が 成り立つ。 K ∂w ∂x x=0 = h0 w|x=0+ h0α− Kβ − h0u¯0 −K ∂w ∂x x=ℓ = hℓw|x=ℓ+ hℓα + (ℓhℓ+ K) β− hℓu¯ℓ (4) よって、α, β を適当に選ぶことで K ∂w ∂x x=0 = h0 w|x=0 −K ∂w ∂x x=ℓ = hℓw|x=ℓ (5) が得られる。これは斉次形の境界条件なので、 変数分離法を適用できる。 例題1 (3) を(2)に代入することで (4)を導 出せよ。2
変数分離法による解法
(3)を満たすw が、熱方程式 wt= kwxx (6) を満たすことは明らかである。uの初期条件を u(x, 0) = φ(x)とすると、wの初期条件は w(x, 0) = φ(x)− α − βx =: φ1(x) と表される。前回と同様の変数分離法より、 w(x, t) = (a cos λx + b sin λx) e−λ2kt (7) が分かる。ここでa, bは任意定数、λは分離定 数である。境界条件(5)より Kλb = h0aKλa sin λℓ− Kλb cos λℓ = hℓa cos λℓ + hℓb sin λℓ
(8) が得られる。よって分離定数は次の方程式を満 たす必要がある。 tan λℓ = K (h0+ hℓ) λ K2λ2− h 0hℓ 端が一定温度にさらされているなら、h0, hℓ → ∞ であることから、a = 0と tan λℓ = 0 すなわち λℓ = nπ が分かる。λ = nπ ℓ を(7)に代入して、重ね合 わせの原理とフーリエ級数展開の原理を利用す ることで w(x, t) = ∞ ∑ n=1 ( 2 ℓ ∫ ℓ 0 φ1(y) sin nπy ℓ dy ) × sinnπx ℓ exp ( −n2π2kt ℓ2 ) (9) が得られる。 演習問題 λ = nπ ℓ を(7)に代入して、重ね合わせの 原理とフーリエ級数展開の原理を利用する ことで(9) を導出せよ（但しa = 0）。 ※ 授業評価アンケートにご協力下さい（9月 15日まで）。うりぼーネット → アンケート
「偏微分方程式」演習問題解答（第
12
回・
2015
年
7
月
13
日）
担当教員：國谷紀良 研究室：工学部本館3W-406 連絡先：tkuniya@port.kobe-u.ac.jp 演習問題 λ = nπ ℓ を(7)に代入して、重ね合わせの 原理とフーリエ級数展開の原理を利用する ことで (9)を導出せよ（但しa = 0）。解答
λ = nπ ℓ を(7) に代入すると（a = 0） w(x, t) = b sinnπx ℓ e −n2π2 ℓ2 kt 重ね合わせの原理より w(x, t) = ∞ ∑ n=1 bnsin nπx ℓ e −n2π2 ℓ2 kt (1) も熱方程式を満たす。t = 0とすると、初期条 件より φ1(x) = ∞ ∑ n=1 bnsin nπx ℓ が得られる。よってフーリエ級数展開の原理 より bn= 2 ℓ ∫ ℓ 0 φ1(y) sin nπy ℓ dy が得られる。これを(1) に代入すれば w(x, t) = ∞ ∑ n=1 ( 2 ℓ ∫ ℓ 0 φ1(y) sin nπy ℓ dy ) × sinnπx ℓ exp ( −n2π2kt ℓ2 ) となり、(9) が得られる。1
数値解法（陽解法）
有限区間x∈ [0, ℓ]の熱方程式 ut= kuxx, k > 0 ux(0, t) = ux(ℓ, t) = 0 (1) を、有限時間 t∈ [0, T ] で考える。J , N を自 然数として、 ∆x := ℓ J, ∆t := T N とし（∆はラプラシアンではない）、 xj := j∆x, j = 0, 1, 2,· · · , J tn:= n∆t, n = 0, 1, 2,· · · , N とすることで、格子点 unj := u (xj, tn) を考えることが出来る。このとき(1)の微分演 算については次の近似が成り立つ。 ut≈ un+1j − unj ∆t , uxx ≈ un j+1− 2unj + unj−1 (∆x)2 また (1)の境界条件は un1 − un0 ∆x = 0, unJ− unJ−1 ∆x = 0 と書き換えられるため、 un0 = un1, unJ = unJ−1 が得られる。以上より、熱方程式(1)は次の離 散的な表現に書き換えられる。un+1j = σunj+1+ (1− 2σ)unj + σunj−1
un0 = un1, unJ = unJ−1 (2) 但し σ := k ∆t (∆x)2 とした。 例題1 上記の手順で熱方程式 (1) を(2)に書 き換えよ。
2
変数分離法による解
変数分離法を利用すると、初期条件u(x, 0) = φ(x)に対する (1)の解は次のようになる。 u(x, t) = ∞ ∑ n=1 ( 2 ℓ ∫ ℓ 0 φ(y) cosnπ ℓ y dy ) × exp ( −kn2π2 ℓ2 t ) cosnπ ℓ x (3)3
数値実験
0 20 40 60 80100 0 2 4 6 8 10-1 0 1 時間 t 位置 x 温度 u(x,t) (a) 0 2 4 6 8 10 -1 -0.5 0 0.5 1 位置 x 温度 u(x,t) t=0 t=25 t=50 t=100 (b) コードの例（MATLAB） k = 1/2; ℓ = 10; T = 100; J = 100; N = 10000; dx = ℓ/J; dt = T/N; s = k∗ dt/(dx2); for x = 1 : 1 : J u(x, 1) = cos(pi∗ x/J); end for t = 1 : 1 : N− 1u(1, t + 1) = s∗ u(2, t) + (1 − s) ∗ u(1, t);
for x = 2 : 1 : J− 1
u(x, t + 1) = s∗ u(x + 1, t)...
+(1− 2 ∗ s) ∗ u(x, t) + s ∗ u(x − 1, t);
end
u(J, t + 1) = (1− s) ∗ u(J, t) + s ∗ u(J − 1, t);
end 演習問題 非斉次項f (x, t) を含む熱方程式 ut= kuxx+ f (x, t), k > 0 を離散的に表現せよ（fjn:= f (xj, tn)）。 ※ 授業評価アンケートにご協力下さい（9月 15日まで）。うりぼーネット → アンケート
「偏微分方程式」演習問題解答（第
13
回・
2015
年
7
月
23
日）
担当教員：國谷紀良 研究室：工学部本館3W-406 連絡先：tkuniya@port.kobe-u.ac.jp 演習問題 非斉次項f (x, t) を含む熱方程式 ut= kuxx+ f (x, t), k > 0 を離散的に表現せよ（fjn:= f (xj, tn)）。解答
un+1j − unj ∆t = k unj+1− 2unj + unj−1 (∆x)2 + f n j となる。これを un+1j について整理すると un+1j = unj + σ(unj+1− 2unj + unj−1) + ∆t fjn= σunj+1+ (1− 2σ)unj + σunj−1+ ∆t fjn
となる。但し σ := k ∆t
