偏微分方程式の線形化
この章では、第1章で扱った弾性体でできた弦の運動方程式を、微小振動の仮定の下に、線形化する。
7.1弦の振動の線形化
この節では、弦の振動が線形の運動方程式を満たす条件を考え、縦波および横波の場合、それぞれ式(2)および式(3)を導く。
7.1.1線形化前
第1章で見たように、3次元空間中に張った弦の運動
𝒙𝑡,𝑥′
について、運動方程式は以下のようになる:(
𝒖
は釣り合い位置からの変位)
𝜌′¨𝒖=𝐸𝛼′(1−𝑃⟂𝛼)𝜕2𝒖𝜕𝑥′2(1)
ただし
𝜌′:釣り合い状態での線質量密度(定数)𝐸:弦のヤング率(定数)𝛼′:釣り合い状態での弦の体積倍率(定数)𝛼𝑡,𝑥′:弦の体積倍率𝑃⟂:弦と垂直な面への正射影行列
である。弦上のラベル
𝑥′
は、釣り合い状態において弦の弧長になっているとする。重力は無視している。
7.1.2線形化:縦波の場合→式(2)
この方程式(1)は線形ではないので、線形化することを考える。もっとも単純なのは、
𝒖
が、常に一直線上に乗っている場合である。この時、
𝑃⟂
がかかる項は消える。よって、以下のような方程式になる:
𝜌′¨𝒖=𝐸𝛼′𝜕2𝒖𝜕𝑥′2
𝒖
のように3次元ベクトルで表す意味があまりないので、
𝒖
の
𝑥
成分座標を、
𝑢∥
で表すことにする:
¨𝑢∥=𝑣2∥𝜕2𝑢∥𝜕𝑥2𝑣∥≡√𝐸𝛼′𝜌′(2)
𝑥=𝑥′
なので置き換えた。これを、(1次元)波動方程式という。
なお、
𝑣
をこのように定義したのは、次節で述べるように、
𝑣
が波動の伝播速度になるからである。
実際に数値計算を行ってみると右図のようになる。初期変位のみを与えた場合である。このように、波動方程式(2)の解は、直感的に予想される通り、振動が伝播していくという波動の性質を示す。特に今の場合、波動の進行方向と変位の方向が平行になっており、これを縦波という。ただし、実際の弾性体がこの方程式を満たすのは、フック弾性体近似が成り立たなければならないので、変位がもっと小さい場合である。
7.1.3線形化:横波の場合→式(3)
縦波以外の場合はどうだろうか。運動方程式(1)を線形化するには、正射影行列
𝑃⟂
が消える必要がある。そのためには、
𝑃⟂
にかかる項が、弦と平行か垂直でなければならない。よって、変位
𝒖
は、弦と平行か垂直である。平行な場合は、上述の縦波の場合なので、垂直な場合を考える。すると、式(1)から
𝑃⟂
が消えて
¨𝒖⟂=𝑣2⟂𝜕2𝒖⟂𝜕𝑥2𝑣⟂≡√𝐸𝛼′𝜌′(1−1𝛼)(3)
となる。線形であるためには、
𝛼
が定数になればよい。これは、弦の伸び縮みが小さいことを意味しているので
弦があまり複雑に波打っていない
という条件を満たせばよい。すると、式(3)は線形になるわけだが、式(2)と同じ形をしている。この場合の波動は、進行方向と垂直なので、横波という。なお、
𝑣⟂<𝑣∥
であるため、振動が伝わる速度は、縦波の場合よりも遅くなる。また、弦を強く張ったほうが、
𝛼
が大きくなるので横波の速度は上がる。
実際に数値計算を行ってみると右図のようになる。初期変位のみの場合である。ただし、上述のように、実際の横波がこの方程式を満たすのは、弦の伸びが無視できるような場合である。
7.21次元波動方程式の初期値問題の解
縦波の方程式(2)と、横波の方程式(3)は、ともに波動方程式
¨𝑢=𝑣2𝜕2𝑢𝜕𝑥2(4)
を満たす。これの初期値問題を解きたい。実際、この方程式は時間の2階微分なので、初期値として初期変位分布
𝑢0,𝑥
と初期速度分布
˙𝑢0,𝑥
を与えれば、その後の変位
𝑢𝑡,𝑥
が決まるはずである。
この節では、波動方程式(4)の初期値問題の解(9)を導く。
7.2.1関数は無限次元ベクトルである
線形偏微分方程式を解くために、すでに扱った線形常微分方程式の議論を使うことを考える。
変数
𝑥
の関数
𝑓(𝑥)
を、
𝑓𝑥
と表しておく。一方、ベクトル
𝒗
の
𝑖
成分は
𝑣𝑖
と表される。両者を見比べると、違いは、「ベクトルの添え字
𝑖
は整数」だが、「関数の添え字
𝑥
は実数」というだけである。即ち、関数は、実数を添え字に持つような、無限次元のベクトルであると見なせる。模式的には、
𝛿𝑥
ごとに刻んで、以下のようになる:
𝑓≈⎡⎢
⎢
⎢
⎢
⎢⎣⋮𝑓−𝛿𝑥𝑓0𝑓𝛿𝑥⋮⎤⎥
⎥
⎥
⎥
⎥⎦,𝒗=⎡⎢
⎢⎣𝑣1𝑣2𝑣3⎤⎥
⎥⎦
𝛿𝑥
を十分小さくしておけば、
𝑓
の関数形をほぼ再現できるので、よい近似になる。
そうすると、微分
𝑑𝑑𝑡
や積分
∫𝑥𝑥′=−∞
などの線形演算子は、無限次元の行列のように見なせそうである:
𝑑𝑑𝑡𝑓≈1𝛿𝑥⎡⎢
⎢
⎢
⎢
⎢⎣⋮𝑓0−𝑓−𝛿𝑥𝑓𝛿𝑥−𝑓0𝑓2𝛿𝑥−𝑓𝛿𝑥⋮⎤⎥
⎥
⎥
⎥
⎥⎦=1𝛿𝑥⎡⎢
⎢
⎢
⎢
⎢⎣⋱⋱−11−11−1⋱⋱⎤⎥
⎥
⎥
⎥
⎥⎦⏟_____⏟_____⏟≈𝑑𝑑𝑥⎡⎢
⎢
⎢
⎢
⎢⎣⋮𝑓−𝛿𝑥𝑓0𝑓𝛿𝑥⋮⎤⎥
⎥
⎥
⎥
⎥⎦∫𝑥𝑥′=−∞𝑓𝑥′≈𝛿𝑥⎡⎢
⎢
⎢
⎢
⎢⎣⋱⋱1⋱11⋱111⋱⋱⋱⋱⋱⎤⎥
⎥
⎥
⎥
⎥⎦⏟____⏟____⏟≈∫𝑥𝑥′=−∞⎡⎢
⎢
⎢
⎢
⎢⎣⋮𝑓−𝛿𝑥𝑓0𝑓𝛿𝑥⋮⎤⎥
⎥
⎥
⎥
⎥⎦
ここで試しに、2つの行列を強引に掛けてみる。すると、ほぼ単位行列になっており、微分と積分がお互いに逆演算であることと整合していることが分かる。とはいえ、微分した後積分すると、積分定数のために元に戻らないので、完全な逆行列にはなってはまずい。これは、無限次元行列の積が非自明であることの現れである。
関数を無限次元ベクトルと見なすように少し見方を変えただけだが、これにより、ベクトルの場合に成り立った性質が、関数に対しても成り立つのではないかという、大胆な類推ができるようになる。このように、関数と線形演算子を、それぞれ無限次元のベクトルと行列とみなして、線形代数のように扱う数学の分野を、関数解析学という。
7.2.21次元波動方程式(5)の初期値問題の解:式(9)
では実際に、この考え方を偏微分方程式(4):(再掲)
¨𝑢=𝑣2𝜕2𝑢𝜕𝑥2(5)
に応用してみよう。
まず、これと同じ形をした常微分方程式は、線形の反発力が働く場合の運動方程式:
¨𝑥=𝛼2𝑥(6)
である。この方程式についてはよくわかっており、初期値問題の解は、力学編第2章で述べたように
𝑥𝑡=(cosh𝛼𝑡)𝑥0+sinh𝛼𝑡𝛼˙𝑥0(7)
となる。
波動方程式(5)と式(6)は似た形をしており、以下の置き換え:
𝛼↔𝑣𝜕𝑥
で一致する。従って、式(5)の解は、解(7)に対して同様の置き換えを形式的に行ったもの:
𝑢𝑡,𝑥=(cosh𝑣𝑡𝜕𝑥)𝑢0,𝑥+sinh𝑣𝑡𝜕𝑥𝑣𝜕𝑥˙𝑢0,𝑥(8)
であると考えられる。とはいえ、赤字部分は、
sinh,cosh
の引数に
𝜕𝑥
が現れたり、
𝜕𝑥
が分母にあったりして見慣れない式である。しかし、これは計算することができるのである。実際、以下の【7.2-注1】のようになる。これを代入すると、波動方程式(5)の初期値問題の解は、初期値を
𝑢0,𝑥,˙𝑢0,𝑥
として、以下のようになる:(ダランベールの解の公式)
𝑢𝑡,𝑥=12(𝑢0,𝑥+𝑣𝑡+𝑢0,𝑥−𝑣𝑡)+12𝑣∫𝑥+𝑣𝑡𝑥′=𝑥−𝑣𝑡˙𝑢0,𝑥′(9)
これが実際に解になっていることを示すのは容易である。波動方程式を満たすことは代入してみれば、各項が波動方程式を満たすことが分かる。初期条件を満たすことは、「
𝑡=0
とすれば右辺の積分は消えて両辺が一致する」こと、および、「両辺を
𝑡
微分してから
𝑡=0
とすると右辺の積分の項だけが残り両辺が一致する」ことから分かる。
なお、解(9)は
𝑢𝑡,𝑥=𝑓−(𝑥+𝑣𝑡)+𝑓+(𝑥−𝑣𝑡)(10)
という形をしている(積分を
∫𝑥+𝑣𝑡𝑥−𝑣𝑡=∫𝑎𝑥−𝑣𝑡+∫𝑥+𝑣𝑡𝑎
と分離してみればよい)。
𝑓−
は、速度
−𝑣
で平行移動する解であり、
𝑓+
は、速度
𝑣
で平行移動する解である。任意の解は、それらの重ね合わせになるわけである。よって、波動の伝播速度が
𝑣
であることが分かる。また、逆に、式(10)の形で表される任意の関数は、波動方程式(5)を満たす(実際に式(5)に代入してみればすぐわかる)。
【7.2-注1】微分演算子 𝜕𝑥 の関数
任意の関数
𝑓𝑥
に、微分演算子
𝜕𝑥
の関数
𝑒𝑎𝜕𝑥,sinh𝑎𝜕𝑥,sinh𝑎𝜕𝑥𝑎𝜕𝑥
を作用させた結果は、それぞれ以下のようになる:
𝑒𝑎𝜕𝑥𝑓𝑥=𝑓𝑥+𝑎(cosh𝑎𝜕𝑥)𝑓𝑥=𝑓𝑥+𝑎+𝑓𝑥−𝑎2sinh𝑎𝜕𝑥𝜕𝑥𝑓𝑥=12∫𝑥+𝑎𝑥′=𝑥−𝑎𝑓𝑥′(11)(12)(13)
指数関数
𝑒𝑎𝜕𝑥
は、
−𝑎
の平行移動を引き起こす演算子なので、
𝜕𝑥
は、平行移動の生成子であるという。
導出
式(11)。以下の微分方程式
˙𝑓=𝑎𝜕𝑥𝑓(14)
の初期値問題の解は
𝑓𝑡,𝑥=𝑓0,𝑥+𝑎𝑡(15)
である(方程式と初期条件を満たすことはすぐ分かる)。式(14)の形式的な解:
𝑓𝑡,𝑥=𝑒𝑎𝑡𝜕𝑥𝑓0,𝑥
(
˙𝒙=𝐴𝒙
の解を
𝒙𝑡=𝑒𝑡𝐴𝒙0
と書くのと同じ)と式(15)の右辺同士を見比べれば、式(11)を得る。
◼
式(12)。
cosh
関数の定義より
cosh𝑎𝜕𝑥=𝑒𝑎𝑡𝜕𝑥+𝑒−𝑎𝑡𝜕𝑥2
である。これに、先ほど求めた式(11)を使えば、よい
◼
式(13)。
sinh𝑎𝜕𝑥
は
(sinh𝑎𝜕𝑥)𝑓𝑥=𝑓𝑥+𝑎−𝑓𝑥−𝑎2
となる。
1𝜕𝑥
は微分の逆演算、即ち、積分である。よって
sinh𝑎𝜕𝑥𝜕𝑥𝑓𝑥=1𝜕𝑥𝑓𝑥+𝑎+𝑓𝑥−𝑎2=∫𝑥+𝑎𝑥′=𝜆𝑓𝑥′−∫𝑥−𝑎𝑥′=𝜆𝑓𝑥′2=12∫𝑥+𝑎𝑥′=𝑥−𝑎𝑓𝑥′
実際には積分定数の自由度が残るが、
𝑎=0
の時に式全体がゼロになるようにとった。
◼
7.3外力がある場合の波動方程式の初期値問題
外力として体積力
𝑓
が働いている場合の波動方程式は
¨𝑢=𝑣2𝜕2𝑢𝜕𝑥2+𝑓(16)
である。この節では、この方程式の初期値問題の解の公式(21)を導出する。
7.3.1偏微分方程式におけるデュアメルの原理
式(16)は、非斉次方程式である。第6章の常微分方程式の場合と同様に、非斉次のデュアメルの原理からこの初期値問題の解を書き下すことができると期待できる。偏微分方程式の場合のデュアメルの原理は、以下の【7.3-注1】で与えられる。
【7.3-注1】偏微分方程式のデュアメルの原理:式(19)
未知の(ベクトル値)関数
𝒖𝑡,𝒙
に対する斉次偏微分方程式:(
̂𝐿
は線形演算子)
˙𝒖𝑡,𝒙=̂𝐿𝒖𝑡,𝒙(17)
において、初期値問題の一般解
𝒖𝑡,𝒙=̂𝐺𝑡𝒖0,𝒙
が求められているとする。
̂𝐺𝑡
は、式(17)の時間発展演算子である。
この時、式(17)にソース項
𝒋𝑡,𝒙
を加えた非斉次方程式:
˙𝒖𝑡,𝒙=̂𝐿𝒖𝑡,𝒙+𝒋𝑡,𝒙(18)
の一般解は、以下のように書ける(デュアメルの原理):
𝒖𝑡,𝒙=̂𝐺𝑡(𝒖0,𝒙+∫𝑡𝑡′=0̂𝐺𝑡′−1𝒋𝑡′,𝒙)(19)
̂𝐺𝑡
を、式(18)における初期値問題のグリーン演算子という。
なお、
̂𝐿
が時間に依存しない場合、
̂𝐺𝑡=𝑒𝑡̂𝐿
となり、式(19)は、以下のように簡単になる:
𝒖𝑡,𝒙=𝑒𝑡̂𝐿𝒖0,𝒙+∫𝑡𝑡′=0𝑒(𝑡−𝑡′)̂𝐿𝒋𝑡′,𝒙(20)
証明
証明は、式(19)を式(18)の左辺に代入するだけである。
◼
第6章の【6.1-注1】で述べたように、常微分方程式の場合にも、同様の公式が成り立つ。
7.3.2外力がある場合の波動方程式の初期値問題の解:式(21)
非斉次波動方程式の場合、デュアメルの原理が使える形(17)にするには、時間微分について1階正規形にする必要がある。
𝑼
および
̂𝐿
を
𝑼≡[𝑢˙𝑢],̂𝐿≡[1𝑣2𝜕2𝑥]
とおけば、1階正規形になる:
𝜕𝑡𝑼=̂𝐿𝑼+[0𝑓]
よって、初期値問題の解は、デュアメルの原理(20)により
𝑼𝑡,𝒙=𝑒𝑡̂𝐿𝑼0,𝑥+∫𝑡𝑡′=0𝑒(𝑡−𝑡′)̂𝐿[0𝑓𝑡′,𝑥]
となる。第1成分だけ分かればよい(第2成分は第1成分を時間微分したもの)。右辺第1項は、斉次の場合の解そのものである。同第2項を計算するには、
𝑒𝑡̂𝐿
の
(1,2)
成分(=右上成分)が必要であるが、これは、式(9)の
˙𝑢0,𝑥
に作用している演算子である。これにより、初期値
𝑢0,𝑥,˙𝑢0,𝑥
およびソース項
𝑓𝑡,𝑥
が与えられている時の初期値問題の解が得られる:
𝑢𝑡,𝑥=12(𝑢0,𝑥+𝑣𝑡+𝑢0,𝑥−𝑣𝑡)+12𝑣∫𝑥+𝑣𝑡𝑥′=𝑥−𝑣𝑡˙𝑢0,𝑥′=+12𝑣∫𝑡𝑡′=0∫𝑥+𝑣(𝑡−𝑡′)𝑥′=𝑥−𝑣(𝑡−𝑡′)𝑓𝑡′,𝑥′(21)
7.3.3シミュレーション
実際に数値計算を行ってみると、初期変位と初期速度をゼロとして、右図のようになる。ただし、実際の弾性体がこの方程式を満たすのは、変位がもっと小さい場合である。