双曲型方程式の例

双曲型方程式の例
1
この講義では、双曲型の偏微分方程式を扱う。最初は、様々な現象に対する波動方程式を導出し、双曲型方程
式の例をいくつか挙げる。
1.1
1.1.1
波動現象の方程式
弦の振動
両端を固定した弦の振動を表す方程式を以下の仮定のもとで導出する。
1. 弦の運動はある平面に限られる(以下それを x 軸と u 軸で定められた平面とする)
2. 弦の各部分は x 軸と垂直な方向にのみ動く(横波)
3. 弦は一様な線密度 ρ を持つ
4. 弦の張力 T は運動の間一定で、弦の接線方向に働く
5. 弦と x 軸とが成す角度 θ は微小である
従って、次の設定である:
独立変数 (independent variable)
:
時刻 t, 1 次元の位置座標 x
未知関数 (unknown function)
:
弦上の点 x の時刻 t での垂直方向の変位 u(t, x)
T
θ(t, x + ∆x)
u(t, x + ∆x)
θ(t, x)
u(t, x)
T
x + ∆x
x
弦の微小部分 [x, x + ∆x] の運動を考え、Newton の第二法則 F = ma を適用する。このとき、微小部分の質量
は m = ρ∆x となり、その加速度は
∂2u
∂t2 (t, x)
となる。また、部分に働く外力 F = (Fh , F ) は、左右両端における
張力の合計である。水平方向の成分 Fh と垂直方向の成分 F は次のように書ける:
Fh
= T cos(θ(t, x + ∆x)) + T cos(θ(t, x))
F
= T sin(θ(t, x + ∆x)) + T sin(θ(t, x))
問題 1.1 cos θ と sin θ を、tan θ を用いて表せ。
tan θ(t, x) = ux (t, x) であるから、
cos(θ(t, x + ∆x)) + cos(θ(t, x)) =
√
√
1/ 1 + (ux (t, x + ∆x))2 − 1/ 1 + (ux (t, x))2
= ∆x ux (t, ξ)uxx (t, ξ){1 + (ux (t, ξ))2 }−3/2
1
がある ξ ∈ (x, x + ∆x) に対し成り立ち、また、
sin(θ(t, x + ∆x)) + sin(θ(t, x)) =
√
√
ux (t, x + ∆x)/ 1 + (ux (t, x + ∆x))2 − ux (t, x)/ 1 + (ux (t, x))2
= ∆x uxx (t, η){1 + (ux (t, η))2 }−3/2
がある η ∈ (x, x + ∆x) に対して成り立つ。
∆x が十分小さければ、Fh ≈ ux F となり、振動の角度が小さいと仮定している (つまり、ux << 1) ので、垂
直方向の成分に比べて水平方向の成分が無視できることがわかる。上の計算より微小部分の運動方程式は次のよ
うになる:
T ∆x uxx (t, η){1 + (ux (t, η))2 }−3/2 = ρ∆x
∂2u
(t, x).
∂t2
この両辺を ∆x で割り、∆x → 0 とすれば、
ρutt (t, x) = T
uxx (t, x)
(1 + ux (t, x))3/2
(1)
という偏微分方程式が得られた。
u も ux も小さいと仮定している。そのため、関数 u を小さいパラメータ ε を用いて、
u(t, x) = εu0 (t, x) + ε2 u1 (t, x) + ε3 u2 (t, x) + . . . ,
のように展開できる。ここで、u0 , u1 , u2 , . . . は t と x について滑らかな関数であるとする。上式を (1) に代入し、
ε の 1 乗の係数をイコールさせると
ρ
∂ 2 u0
∂ 2 u0
(t,
x)
=
T
(t, x)
∂t2
∂x2
(2)
が得られる。
問題 1.2 ε の 2 乗と 3 乗の係数より従う微分方程式を書け。
ε を十分小さくして、ε2 u1 + ε3 u2 + . . . が εu0 に比べて無視できる程度にすれば、弦の運動を εu0 のみで記述
できることになる。よって、u(t, x) = εu0 (t, x) とすると (2) より u が次の方程式を満たすことがわかる:
ρ
定数 c =
√
T /ρ を導入して、上式は
∂2u
∂2u
(t,
x)
=
T
(t, x).
∂t2
∂x2
utt − c2 uxx = 0 の形に書かれ、空間 1 次元の波動方程式 (wave equation) と呼ばれる。
注意: 導出を簡単にするため多くの仮定をおいて得られた方程式なので、この「モデル」が現象を十分に表
現できているかを数学的に・物理的に確認しなければならない。
外力 (outer force)
以上で考えた弦に重力が作用するとすれば、F に −gρ∆x が加わる。より一般に、時
刻 t において弦上の点 x で弦に垂直な方向に単位質量当り f (t, x) の力が働くとき、f (t, x)ρ∆x が加わり、得られ
る偏微分方程式は
∂2u
∂2u
= c2 2 + f
2
∂t
∂x
となる。
抵抗 (resistance)
弦が空気から抵抗を受けることを考えるとき、最も基本的な抵抗は速度に比例するも
のである。速度抵抗を取り入れた方程式は次のようになる。
2
∂u
∂2u
2∂ u
+
k
−
c
=0
∂t2
∂t
∂x2
2
1.1.2
バネの振動
波は基本的には 2 種類ある:
• 横波: 波の進行方向と垂直方向に媒質が振動する波(例えば、弦の振動)
• 縦波: 波の進行方向と同一方向に媒質が振動する波(例えば、バネの振動)
この節では、縦波の例であるバネの振動を考える。
力がかからないときのバネの長さを l とし、バネが伸ばされ区間 [a, b] の端点に対応する点 A と B でそれぞ
れの先端を固定されているとする(下図を参照)。
b
A
b
P
( )
u t; x
b
a
x
b
( + )
u t; x
B
b
x
+ x
b
x
Q
b
長さ d のバネにかかる外力とバネの伸長の関係がフックの法則により与えられる。つまり、外力は伸びに比例
する:
s
F =k .
d
今回もバネの微小部分 ∆S = [x, x + ∆x] を考える。この両端が点 P と Q に対応するとする(図を参照)。また、
α = l/(b − a) とおく。バネが伸ばされ固定される前の自然な状態のときの微小部分の長さは ∆x · l/(b − a) = α∆x
となる。
今、時間 t が経ったあと、点 x が u(t, x) だけ移動するとする。すなわち、点 P は座標 x + u(t, x) に移り、点
Q は座標 x + ∆x + u(t, x + ∆x) に移動する。よって、時刻 t において微小部分 ∆S の長さは
∆x + u(t, x + ∆x) − u(t, x)
となり、固定されていない状態と比べて
u(t, x + ∆x) − u(t, x) + (1 − α)∆x
だけ伸ばされたことになる。
フックの法則より ∆S の両端において発生する外向きの力が
k
u(t, x + ∆x) − u(t, x) + (1 − α)∆x
α∆x
になることが従う。∆x → 0 のとき、上の比が
{
(1
)}
1 ∂u
(t, x) +
−1
k
α ∂x
α
に収束し、この式は P における力を表している。∆x > 0 に対し、∆S の両端に働く力を書くと、
k(ux (t, x)/α + 1/α − 1)
k(ux (t, x + ∆x)/α + 1/α − 1)
...
...
左方向に働く P における力
右方向に働く Q における力
となる。よって、∆S 全体には右方向の力
k
ux (t, x + ∆x) − ux (t, x)
k
= ∆x uxx (t, ξ)
α
α
3
が働く。ただし、ξ ∈ (x, x + ∆x) である。
区間 AB に伸ばされたバネの密度を ρ とおくと、ニュートンの第 2 法則より
ρ∆x utt (t, x) =
k
∆x uxx (t, ξ)
α
が得られる。両辺を ∆x で割り、∆x → 0 とすれば、
utt (t, x) =
k
uxx (t, x)
αρ
が従う。主な動的要因が弦の運動とバネの運動とで異なるにもかかわらず、最終的な運動方程式が同じ形になる
ことが興味深い事実である。
以下の説明で様々なタイプの作用素が出てくるので、ベクトル解析の基本的な知識を復習しておく。
復習.
基礎的な演算子 ”div”, ”∇”, ”rot”と”∆”.
• a = (a1 , . . . , an ) というベクトル場の ”div” または ”∇·”(発散)は、流れに
よって生じるある点での流量の発生(体積変化)を表すスカラーである:
div a =
n
∑
∂ai
∂xi
i=1
(
div a = div (a1 , a2 ) =
∂a1
∂a2
+
∂x1
∂x2
)
2 次元の場合
• スカラー関数 u の ”∇” または ”grad”(勾配)は、大きさはその点での最大
の傾斜の大きさで、方向は等高線に垂直であるようなベクトルである:
∇u =
(
( ∂u ∂u
∂u )T
,
,...,
∂x1 ∂x2
∂xn
∇u =
( ∂u ∂u )T
,
∂x1 ∂x2
)
2 次元の場合
• ベクトル場 v の”rot” または ”∇×”(回転)は大きさは面積当りの渦の量で、
方向は渦がある面と垂直にとるようなベクトルである(普通は 3 次元の場合だ
け考える)
:
e1
e2
e3 ( ∂v3 ∂v2 ∂v1 ∂v3 ∂v2 ∂v1 )T ∂
∂
∂ rot v = rot (v1 , v2 , v3 ) =
−
,
−
,
−
= ∂x1 ∂x
∂x3 2
∂x2 ∂x3 ∂x3 ∂x1 ∂x1 ∂x2
v1
v2
v3 • ”rot”作用素は”curl”と呼ばれることもあるが、ここでは”curl”という記号を 2
変数のスカラー関数に限定して使い、次のように定義する
curl v = (
∂v
∂v T
,−
)
∂x2
∂x1
• ”∆”はラプラシアンと言われる作用素で、スカラー関数 u に対し次のように定
義される:
∆u =
n
∑
∂2u
i=1
∂x2i
(
∆u =
4
∂2u ∂2u
+ 2
∂x21
∂x2
)
2 次元の場合
例.
関数 f が 2 回微分可能な 3 変数のスカラー関数であるとする。そのとき rot (∇f ) は零ベクトルである。
証明.
[( ∂f ∂f ∂f )T ]
,
,
∂x1 ∂x2 ∂x3
( ∂ ( ∂f )
∂ ( ∂f ) ∂ ( ∂f )
∂ ( ∂f ) ∂ ( ∂f )
∂ ( ∂f ))T
=
−
,
−
,
−
∂x2 ∂x3
∂x3 ∂x2 ∂x3 ∂x1
∂x1 ∂x3 ∂x1 ∂x2
∂x2 ∂x1
( ∂2f
)
2
2
2
2
2
T
∂ f
∂ f
∂ f
∂ f
∂ f
=
−
,
−
,
−
∂x2 ∂x3
∂x3 ∂x2 ∂x3 ∂x1
∂x1 ∂x3 ∂x1 ∂x2
∂x2 ∂x1
rot (∇f ) =
rot
関数 f の 2 階微分が連続なので、微分する順序を好きに変えてもよく、結果が直ちに従う。
問題 1.3 上の作用素に関する説明を参考にして、以下の問いに答えよ。
(a) 流体の速度場が式 v(x, y, z) = (x sin2 y, z − y, 1 + z cos2 y) で与えられるとき、流体が非圧縮かどうか判断
せよ。
(b) 関数 f (x, y) = x2 e1−y が点 (x0 , y0 ) = (2, 1) において最も急に減少する方向を求めよ。
(c) ベクトル場 v(x, y, z) = (−y, x, 0) に対し原点における回転ベクトルを計算せよ。
(d) 関数 f (x, y) = 1 + x2 − y 2 が調和関数であるかどうか判断せよ(関数 f のラプラシアン ∆f が領域のすべ
ての点で零になるとき、 f がその領域で調和関数であるという)。
問題 1.4 3 次元空間の場合、次の関係式が成り立つことを確認せよ。ここで、f はスカラー関数、u, v は 3 次元
のベクトル関数である。(”∇·”は”div”と同じ意味で、”∇×”は”rot”のことをさす)
(A)
∇ · (f v) = f ∇ · v + v · ∇f
(B)
∇ · (u × v) = v · (∇ × u) − u · (∇ × v)
(C)
∇ × (f u) = ∇f × u + f ∇ × u
(D)
∇ · (∇ × u) = 0
(E)
∇ × (∇ × u) = ∇(∇ · u) − ∆u
5
ベクトル解析の基本定理も使うので、ここで復習しておく。
復習.
グリーンの定理について
ガウスの発散定理
∫∫∫
∫∫
u · n dS
div u dV =
R
(3)
S
物理的な意味:ある閉領域 R 内での全発散は閉領域の表面から外へ流れ出る全流量
に等しい。
u
n
rot v
n
rot v · n
u·n
dS
D
R
v·t
v
∂D
ガウスの定理
ストークスの定理
ストークスの定理
∫∫
∫
(rot v) · n dS =
D
v · ds
(4)
∂D
物理的な意味:閉領域 D 内の全循環は閉領域の周囲の閉曲線 ∂D 上の線積分に等し
い。
グリーンの定理
∫∫
∫∫
I
ϕ∆ψ dx dy = −
Ω
∂ψ
ds
∂n
(5)
ϕ(a · n) ds
(6)
∇ψ · ∇ϕ dx dy +
Ω
ϕ
∂Ω
(ガウスの定理で u = ϕ∇ψ とおくと得られる)
さらに、ガウスの定理で u = ϕa とすると、
∫∫
∫∫
I
div a ϕ dx dy +
a · ∇ϕ dx dy =
Ω
Ω
∂Ω
という形のグリーンの定理が得られる。
問題 1.5 ベクトル場 f を次のように定義する:
f (x, y, z) = (y, x − 2xz, −xy).
このベクトル場に対し、
∫
(rot f ) · n dS S
の値を計算せよ。
ここで、S は球面 x2 + y 2 + z 2 = r2 , r > 0 で、ベクトル n は S への外向き単位法線ベクトルである。
6
1.1.3
音の伝播
本節では、音の伝達を表す方程式を導出する。簡単のため、外力が働かない場合を考える。
1. ニュートンの第 2 法則: 気体の任意の部分領域 D に対して述べる。ニュートンの法則に登場する力は、ま
わりの気体が領域 D に及ぼす圧力になり、次のように書ける:
∫
− p(x)ν(x) dS
S
ここで、S = ∂D、p(x) は気体の点 x における圧力、ν ∈ R3 は外向き単位法線ベクトル、そして dS は S
の面積要素である。発散定理を用いると、上の積分は以下の形に直される:
∫
−
∇p(x) dx.
D
時刻 t と点 x における気体の密度を ρ(t, x) とし、気体の加速度を w(t, x) とする。そのとき、ニュートン
∫
の第 2 法則は
∫
−
∇p(x) dx =
ρ(t, x)w(t, x) dx
D
D
という形をとる。この式はすべての領域 D に対して成り立つから、以下の式が任意の t と x について成り
立つことがわかる:
ρ(t, x)w(t, x) = −∇p(t, x).
(7)
2. 加速度の式: ある点と時刻における気体の速度を v(t, x) という記号で表す。時刻 t において点 x にあった気
体の粒子は時間 ∆t 後に点 x + ∆tv(t, x) に移る。従って、時間 ∆t 後の粒子の速度は v(t + ∆t, x + ∆tv(t, x))
となる。
問題 1.6 加速度 w を速度 v の式で表せ。
よって、
w(t, x) =
∂v
(t, x) + ∇v(t, x) · v(t, x).
∂t
(8)
3. 連続の式: 次に速度 v と密度 ρ の関係を求めたい。この関係は気体の任意の部分領域 D に対し質量保存
∫
の式を書くことにより明らかにされる。領域 D における気体の質量が D ρ(t, x) dx であるので、D のなか
∫
の質量の変化率は
D
∂ρ
(t, x) dx.
∂t
となる。この量は境界 S を通して単位時間当たりに D から出る気体の質量に打ち消されるから、
∫
∫
∫
(
)
∂ρ
(t, x) dx = − ρ(t, x)v(t, x) · ν(x) dS = −
div ρ(t, x)v(t, x) dx.
∂t
S
D
D
上の変形で再び発散定理を用いた。
よって、
∫
D
∂ρ
(t, x) dx = −
∂t
∫
(
)
div ρ(t, x)v(t, x) dx.
D
上式が任意の領域 D について成立するので、
(
)
∂ρ
(t, x) + div ρ(t, x)v(t, x) = 0.
∂t
この式は普段、連続の式と呼ばれる。
7
(9)
上の準備を利用して音の伝播の方程式は容易に導かれる。定常状態の空気の密度を ρ0 とし、空気の圧縮率を
表す関数として u(t, x) を導入する。すなわち、
ρ(t, x) = ρ0 (1 + u(t, x)).
今、圧力が密度に比例すると仮定する。つまり、体積膨張率とよばれる定数 K が存在し、
p(t, x) = Kρ(t, x) = Kρ0 (1 + u(t, x))
が成り立つ。この関係式より次を得る:
∇p = Kρ0 ∇u.
(10)
空気がわずかにしか動かない状況を考えたいので、小さいパラメータ ε > 0 を用いて次のようにおく:
v
= εv 0 + ε2 v 1 + . . . ,
u
= εu0 + ε2 u1 + . . . .
式 (8) と式 (10) を質量保存の式 (7) に代入し、さらに v と u を展開した形を代入して、ε の同じベキを比較す
る。ε の 1 乗に対し、次の関係式が得られる
∂v 0
= −Kρ0 ∇u0 .
∂t
(11)
∂u0
+ ρ0 divv 0 = 0.
∂t
(12)
ρ0
同様の方法で連続の式 (9) より次を得る:
ρ0
問題 1.7 式 (11) で発散をとり、式 (12) で時間についての偏微分をとり、得られる方程式から速度 v 0 を消去
せよ。
計算すると、以下の関係式にたどり着く:
∂ 2 u0
− K∆u0 = 0.
∂t2
この方程式を解けば、v 0 は式 (11) を用いて計算できる。さらに、高位近似 uj , v j (j = 1, 2, . . . ) も順番に求める
ことができる。
上の考察より、振幅の小さい音波の伝播を考えるとき、ε2 以上の項が無視してもよいことになる。従って、u
は方程式
∂2u
− K∆u = 0
∂t2
√
に支配される。K が音波の速度に関係があり、すなわち K が電波速度であるということが後で分かる。
問題 1.8 体積膨張率 K が定数ではなく空間の位置に依存するときの音の伝播の方程式を導出せよ。
K が一定でないとき、音波は直線上ではなく K(x) に依存する曲線上に伝播する。そのため、例えば、地球の
表面が冷え空気が暖かいときに遠くの音が聞こえることがある。
8
1.1.4
マクスウェル方程式
マクスウェルの方程式は電荷と電流がどのように電場と磁場を生成するかについて記述する連立方程式である。
• ガウスの式と磁束保存の式は電荷から発生する電場・磁場を表している。
• アンペールの式は電流や変動電場の周りの磁場の振る舞いを表している。
• ファラデーの式は変動磁場の周りの電場の振る舞いを表している。
ε0 div E
rot E
div B
1
rot B
µ0
= q,
∂B
,
= −
∂t
= 0,
∂E
= ε0
+ j,
∂t
(ガウスの式)
(13)
(ファラデーの式)
(14)
(磁束保存の式)
(15)
(アンペールの式)
(16)
記号の意味:
E
...
電場の強度 [V/m]
B
j
q
...
...
...
磁束密度 [Wb/m2 ]
電流密度 [A/m2 ]
電荷密度 [C/m3 ]
µ0
...
透磁率
ε0
...
誘電率
q と j が与えられたとして、電場 E のそれぞれの成分と磁場の B のそれぞれの成分は計 6 個の未知関数に
なる。
次の性質を満たすローレンツゲージという組 (A, ϕ) が存在することが知られている。つまり、A は R3 値関数
で、ϕ は実数値関数で、次の関係式を満たす:
ε0
B
= rot A
(17)
E
∂A
= −∇ϕ −
∂t
(18)
∂ϕ
1
+
div A =
∂t
µ0
0.
(19)
このことを使って、マクスウェル方程式が事実上、波動方程式であることを示す。
問題 1.9 式 (18) をガウスの式に代入し、その結果を式 (19) と組み合わせることにより ϕ に対する以下の方程
式を導け:
ε0 µ0
∂2ϕ
q
− ∆ϕ = .
∂t2
ε0
(20)
問題 1.10 式 (17) と式 (18) をアンペールの式に代入し、その結果を式 (19) の勾配と組み合わせることにより
A の以下の方程式を導け:
∂2A
− ∆A = µ0 j.
∂t2
計算では ∇div − rot rot = ∆ が成り立つことを用いる。
ε0 µ0
(21)
逆に、(A, ϕ) が (20) と (21) を満たし、関係式 (19) も成立するとき、(17) と (18) それぞれにより与えられる
B と E がマクスウェル方程式を満たすことが簡単に示される。
9
1.1.5
弾性方程式
ここでは弾性体の中で波が伝播する現象を表す方程式を紹介する。例えば、地球の地殻は弾性体であり、地震
はそのなかを伝わる弾性波にすぎない。
Ω ⊂ R3 が定常状態にあり力の働かない弾性体の一部であるとする。また、時刻 t において点 x が u(t, x) だ
け移動するとする。この変位を表すベクトル関数 u(t, x) が次の弾性方程式と呼ばれる式を満たすことが知られ
ている:
∂2u
= α∆u + β∇div u + f .
∂t2
ここで、f は外力、ρ は弾性体の密度である。なお、定数 α = µ, β = µ + λ は弾性体の性質を表すラメ定数 λ, µ
ρ
により定められる。
一様な地質がある場所で地震が起こるとき、地震波の伝播は弾性方程式により記述される。そのような場合、
違う性質を持つ 2 種類の波が現れる:
√
(α + β)/ρ の縦波
√
• 電波速度 α/ρ の横波
• 伝播速度
縦波の伝播速度が横波の速度より大きいので、縦波がある地点に先に伝わり、P-波 (primary wave) と呼ばれる。
一方、後からくる横波はふつう S-波 (secondary wave) と呼ばれる。
1.2
保存則系
応用でよく現れる双曲型方程式のもうひとつのグループは 1 階微分方程式で、「保存則」と呼ばれるものであ
る。1 次元の保存則は発散型の特別な非線形微分方程式系で、
∂u
∂
+
F (u) = 0
∂t
∂x
の形をとる。ここで、(x, t) ⊂ R2 は平面上の座標、u : R2 → RN は未知のベクトル値関数、そして F : RN → RN
は与えられた滑らかな関数である。
u ∈ C 1 であれば、この系は以下の同値の形に書き換えられる:
∂u
∂u
+ F 0 (u)
= 0.
∂t
∂x
ここで、行列 F 0 (u) の固有値がすべて実数で互いに異なると仮定される。
最も簡単で有益な例はバーガース方程式である:
∂u
∂ ( u2 )
+
= 0.
∂t
∂x 2
また、物理的な例としてオイラー方程式系が最も有名である:
ρt + (ρu)x
2
(
(ρ(u2 /2 + e)
(ρu)t + (ρu + p)x
)
+ ρu(u2 /2 + e + p/ρ)
)
(
t
=
0
=
0
=
0.
x
ここで、ρ は流体の密度、u は流体の速度、p は圧力、そして e は内部エネルギーを意味する。関数 e = e(ρ, p)
は流体の性質についての物理的な考察により与えられるため、上の 3 本の式において (ρ, u, p) が未知関数となる。
10
参考文献
1. Mitsuru Ikawa, Hyperbolic Partial Differential Equations and Wave Phenomena, AMS, 2000. (日本語
版もある)
2. Peter D. Lax, Hyperbolic Partial Differential Equations, AMS, 2006.
3. Serge Alinhac, Hyperbolic Partial Differential Equations, Springer, 2009.
4. Lars H¨ormander, Lectures on Nonlinear Partial Differential Equations, Springer, 1997.
5. Peter D. Lax, Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves,
SIAM, 1973.
6. Lawrence C. Evans, Partial Differential Equations, AMS, 2008.
7. Alfio Quarteroni, Alberto Valli, Numerical Approximation of Partial Differential Equations, Springer,
2008.
11