第 9 章 最後に応用を 2 つ 110
F.4 Z 上の関数の離散時間 Fourier と畳み込み
G.1.2 波動方程式の初期値問題 , d’Alembert の波動公式
(G.1) を満たす解 u=u(x, t) は無数に存在し、解を一つに特定するためには、何か他に条件を追
加しなければならない。代表的なものは、特定の時刻 (ここでは t= 0) における状態を指定する条 件である。
(G.4) u(x,0) =ϕ(x), ∂u
∂t(x,0) =ψ(x) (x∈R)
条件 (G.1) と (G.4) を同時に満たす u = u(x, t) を求めよ、という問題を初期値問題 (initial value problem)またはCauchy 問題(Cauchy problem)と呼ぶ。また(G.4)を初期条件(initial condition)と呼ぶ。波動方程式は時刻 tについて 2階であることから、初期条件の式が二つになっ ていることに注目しよう。二つの関数 ϕ と ψ を初期値(initial values)と呼ぶ。
前項で見たように、(G.1) を満たす uの一般形 (G.2) が得られた。そこで、初期条件 (G.4)を満 たすように、(G.2) 中の f, g を定められれば初期値問題が解けることになるが、これが実際に可能 であることを以下に示す。初期条件に (G.2)を代入すると
f(x) +g(x) = ϕ(x) (x∈R), c(−f′(x) +g′(x)) =ψ(x) (x∈R) となる。後者から
−f(x) +g(x) = 1 c
∫ x
0
ψ(y)dy−f(0) +g(0) (x∈R) が得られるので、連立1次方程式を解いて
f(x) = 1 2
(
ϕ(x)− 1 c
∫ x
0
ψ(y)dy+f(0)−g(0) )
, g(x) = 1
2 (
ϕ(x) + 1 c
∫ x
0
ψ(y)dy−f(0) +g(0) )
. これから
(G.5) u(x, t) = 1
2(ϕ(x−ct) +ϕ(x+ct)) + 1 2c
∫ x+ct
x−ct
ψ(y)dy.
これを ダ ラ ン ベ ー ル
d’Alembert の波動公式あるいはStokes の公式と呼ぶ (Euler による、1750年頃の結果4)。
以上をまとめると次の定理を得る。
定理 G.1.2 (1次元波動方程式の初期値問題の解の一意存在) ϕ ∈ C2(R), ψ ∈ C1(R) ならば、
波動方程式の初期値問題 (G.1), (G.4)の C2 級の解が一意的に存在し、それは d’Alembert の波
動公式 (G.5) で与えられる。
4Eulerが得た結果なのに、なぜd’AlembertやStokesの名を冠されるのだろうか？ともあれ、簡単であるから1750 年という早い時期に解けたのだろうが、偏微分方程式の問題がこんなに簡単かつ鮮やかに解けるのは奇跡のような気が する。
興味があったらやってみよう 次の各場合に d’Alembert の波動公式で解を計算し、それがどうい う波を表わすか考えてみよ (コンピューターが使える場合、アニメーションを表示してみよう)5。 (1) ϕ(x) = sinx, ψ(x)≡0.
(2) ϕ(x) = {
1 (x∈(−1,1))
0 (それ以外) , ψ(x)≡0.
(3) ϕ(x) = 0, ψ(x) = sinx.
問 9. ϕ ∈C2(R;R), ψ ∈C1(R;R) とするとき、
u(x, t) := 1
2(ϕ(x−ct) +ϕ(x+ct)) + 1 2c
∫ x+ct
x−ct
ψ(y)dy ((x, t)∈R×R) で定めた u が
1
c2utt(x, t) =uxx(x, t) in R×R, u(x,0) =ϕ(x) (x∈R), ut(x,0) =ψ(x) (x∈R) を満すことを直接計算で示せ。
5実は(2)の場合、d’Alembertの波動公式に機械的に代入することによってuは簡単に「求まる」が、不連続関数に なってしまう。不連続関数をどうやって微分方程式の解と考えるのか、頭をひねらねばならないところであるが、この 結果は物理的にそれなりにもっともなものになっている。この辺りの事情は後述する熱伝導方程式とは対照的である。
熱伝導方程式の解は、最初に不連続性があっても、すぐに解は滑らかになってしまうが、波動方程式の解においては、
不連続性は消滅せずに伝搬していく。実用面から見ても、波動現象においては、衝撃波のような、不連続な関数で表現 するほうがむしろ適当なものを扱う必要があるため、これは重要な問題である。
付 録 H [0, ∞ ) 上の関数の Laplace 変換と 畳み込み
f: [0,∞)→C とするとき、fe:R→C を f(x) :=e
{
f(x) (x≥0) 0 (x <0)
で定めて、f の代わりに feを用いて、Fourier変換や畳み込みを定義することも出来るが、Laplace 変換を使うという方法も有用である。
f: [0,∞)→C とするとき、
Lf(s) :=
∫ ∞
0
f(x)e−sxdx (s は複素変数) で定まる Lf を f の Laplace 変換(the Laplace transform of f) と呼ぶ。
g の 逆 Laplace 変換(the inverse Laplace transform ofg) とは Lf =g
を満たす関数f のことをいう。存在するならば一意に定まることが知られている。f =L−1g と書 くことにする。
余談 H.0.1 (Laplace 逆変換の積分表示) g の特異点を含まない領域{z ∈C|Rez > γ}において、
g の Laplace 逆変換L−1g は、Bromwich 積分 L−1g(t) = lim
R→+∞
1 2πi
∫ γ+iR
γ−iR
g(s)estds で求められる。
f, g: [0,∞)→C のとき、f∗g は f ∗g(x) =
∫ x 0
f(x−y)g(y)dy (x∈[0,∞)) で定義される関数 f∗g: [0,∞)→C である。実は
f∗g(x) =fe∗eg(x) (x≥0)
であるから、?? 節の結果から、こちらの畳み込みの性質のいくつかを導くことが出来る。
f と g の畳み込みのLaplace 変換は、f と g の Laplace 変換の積である: L[f ∗g] =LfLg.
付 録 I Fourier 級数の部分和のグラフ
I.1 問
(h については、デルタ関数を知っている人だけ解答せよ。) 次の 3つの関数の Fourier 級数を求 め、コンピューターを用いて部分和のグラフを描け(何項取るかは、いくつか試してから自分で決め て)。f, g,h の関係について気づいたことがあれば説明せよ。
f(x) =|x| (x∈[−π, π]), (f は周期2π), g(x) =
1 (x∈(0, π)), 0 (x= 0,±π),
−1 (x∈(−π,0))
(g は周期 2π), h(x) = 2∑
n∈Z
(−1)nδ(x−nπ) (δ は Diracのデルタ関数).
関数 f, g の Fourier 級数展開そのものは授業中に説明した通り、
f(x) = π 2 − 4
π
(cosx
12 +cos 3x
32 + cos 5x 52 +· · ·
) , g(x) = 4
π
(sinx
1 +sin 3x
3 + sin 5x 5 +· · ·
) .
f は連続で区分的にC1級、g は区分的に C1 級 (ゆえに2乗可積分) であるから、どのような収束 をするかは、10/17 の授業で説明した3つの定理で説明が出来る。
こうすれば Fourier 級数の部分和fn, gn のグラフが描ける
f[n_,x_]:=Pi/2-4/Pi Sum[Cos[k x]/k^2,{k,1,n,2}]
gf10=Plot[f[10, x], {x, -10, 10}]
g[n_,x_]:=4/Pi Sum[Sin[k x]/k,{k,1,n,2}]
gg20=Plot[g[20,x],{x,-10,10}]
nを変えてgn の変化の様子を見たければ、
Manipulate[Plot[g[n, x], {x, -Pi, Pi}, PlotPoints -> 200], {n, 1, 100, 1}]
グラフィックスを保存するには Export["graphf10.pdf",gf10]
Export["graphg20.pdf",gg20]
g の不連続点 x = 0,±π の近くでは、ぴょこっと高い山、それと対称な低い谷があり、n が増加 すると横幅は小さくなるが、高さ(深さ)はほとんど変わらない (小さくなってくれない) ことが見 て取れる。
おまけ: hn のグラフを描いてみる
h[n_,x_]:=4/Pi Sum[Cos[k x],{k,1,n,2}]
Manipulate[Plot[h[n,x], {x,-7,7}, PlotRange -> Full, AspectRatio->Automatic, PlotPoints -> 100], {n, 2, 20, 2}]
, AspectRatio->Automatic は入れた方が良いかどうか判らないが。
もっと色々 Mathematica にさせてみる
f やg のグラフを描かせたり、Fourier係数を求めて、部分和の計算をさせたりしてみよう。f(x) =
|x| (−π < x < π) ということで、積分の計算をするだけならば
f[x_]:=Abs[x]
使うだけで、十分であるが、周期関数らしいグラフを描くには次のようにすると良い。
f[x_]:=Abs[Mod[x,2Pi,-Pi]]
gf=Plot[f[x],{x,-10,10}]
Fourier 係数ak = π1 ∫π
−πf(x) coskx dx の計算は、例えば
Integrate[f[x]Cos[k x],{x,-Pi,Pi}]/Pi
これで2(−1+coskπ+kπsinkπ)
k2π が得られる。k が整数と教えてないので、結果が複雑。
Simplify[Integrate[f[x]Cos[k x],{x,-Pi,Pi}]/Pi, Assumptions->Element[k,Integers]]
とすると、2(−1+(−1)k)
k2π が得られる。実はこれにも問題がある。k= 0 のとき正しくない式である！！ ここら辺は Mathmatica は変な割り切り方をしている(マニュアルによると、“不定積分において，
Integrateは，ほとんどすべてのパラメータの値に対して正しい結果を求める．” —少数の例外は無
視？)。注意が必要である。k = 0 の場合については別扱いで
Integrate[f[x],{x,-Pi,Pi}]/Pi,
とすると良い。π という結果が得られる。
bk= 1π∫π
−πf(x) sinkx dxについては
Integrate[f[x] Sin[k x], {x, -Pi, Pi}]
とすると 0が返る。
f のFourier 係数をa[] という関数で計算させて、f の n 項までの部分和 fn を計算するように
する。
ほぼ全自動
Clear[a,f]
f[x_]:=Abs[Mod[x,2Pi,-Pi]]
a[k_]=Simplify[Integrate[f[x]Cos[k x],{x,-Pi,Pi}]/Pi, Assumptions->Element[k,Integers]]
a[0]=Integrate[f[x],{x,-Pi,Pi}]/Pi
f[n_,x_]:=a[0]/2+Sum[a[k]Cos[k x],{k,1,n}]
gf10=Plot[f[10,x],{x,-10,10}]
(実は Simplify[] はしなくてもあまり効率に影響がない。f[n,x] は Sum[]を使っているので、遅
延評価 := をしないとうまく動かない。一方 a[k ]の方を遅延評価するのはお勧めできない。ここ ら辺は難しい。最初に Fourier 係数を全部計算して、その結果をリストに保存して使うのが良いか もしれない。)
x= 0 のところでfn(x)−f(x)の値を調べると、図 I.1のように減少している(1/n に比例してい ることが分かる)。
list=Table[f[2^k,0.0]-f[0.0],{k,1,10}]
gerror=ListLogPlot[list]
Export["error.pdf",gerror]
n = 210 のとき誤差は 0.000621699≒0.6×10−3
2 4 6 8 10
10-3 10-2 0.05 0.10 0.50 1
図 I.1: fn(0)−f(0) の変化(n = 2k,k = 1,2,· · · ,10)
Clear[g,bg]
g[x_]=Sign[Mod[x,2Pi,-Pi]]
bg[k_]=Simplify[Integrate[g[x]Sin[k x],{x,-Pi,Pi}]/Pi, Assumptions->Element[k,Integers]]
g[n_,x_]:=Sum[bg[k]Sin[k x],{k,1,n}]
ggandg10=Plot[{g[x],g[10,x]},{x,-10,10}]
-10 -5 5 10
-1.0 -0.5 0.5 1.0
図 I.2: g と g の10項までの部分和 g10 のグラフ
h について
h(x) = 2∑
n∈Z
(−1)nδ(x−nπ).
ここで δ は Dirac のデルタ関数 (以下では単にデルタ関数という) を表す。デルタ関数は実は関
数ではなく、超関数というものであり、本来はグラフを描いたりすることも出来ないものであるが (x̸= 0 で 0, x= 0 で ∞ という「説明」もあるが、∞ は図に描けない)、実は Fourier 級数展開は 可能であり、その部分和は以下に見るように普通の関数なので、グラフを描くことが出来る。
デルタ関数は、任意の連続関数 φに対して
∫ ∞
−∞
δ(x)φ(x)dx=φ(0) を満たすという条件や、
δ(x) = 0 (x̸= 0),
∫ ∞
−∞
δ(x)dx= 1
で特徴づけられる。物理学では、大きさのない1点に質量や電荷などが集中しているという理想化 を考えるが(質量の場合は質点、電荷の場合は点電荷と呼ぶ)、デルタ関数は、そのように1点に1 単位のものが集中している場合の、密度を表す関数と考えられる(直観的には、原点で無限大。積分 すると総量になって 1)。
h の定義には、デルタ関数を平行移動した δ(x−nπ) が現れるが、δ(x−a) は
∫ ∞
−∞
δ(x−a)φ(x)dx=φ(a) という性質を持つ。
さて、h の Fourier 級数展開の計算であるが、f や g のように [−π, π] で積分するとして、どう
すれば良いだろうか？
(♯1) h(x) = 2 (δ(x)−δ(x−π)),
(♯2) h(x) = 2 (δ(x)−δ(x+π)),
(♯3) h(x) = 2 (δ(x)−δ(x−π)−δ(x+π)),
このどれが正しいのだろうか？周期2πの周期関数とは R/2πZで定義された関数のことと考えると、
π と −π は同じ点を表すことになるので、(♯1), (♯2)のどちらも同じことを表していて正しそうであ り、(♯3) は1つのものを重複してカウントしていて間違っていそうである。
あるいは積分範囲を [−π, π]でなく、少し右にずらした [−π+ε, π+ε]を考えてみると (εは小さ い正の数)、(♯1)を採用することになり、逆に [−π−ε, π−ε]と左にずらすと、(♯2) を採用すること になる。どちらで計算しても Fourier 係数は同じになる。
これから
ak = 1 π
∫ π
−π
h(x) coskx dx= 2
π(cos (k·0)−coskπ) = 2 π
(1−(−1)k) , bk= 1
π
∫ π
−π
h(x) sinkx dx= 2
π(sin (k·0)−sinkπ) = 0.
ゆえに h の Fourier級数展開は h(x) = 4
π(cosx+ cos 3x+ cos 5x+· · ·).
この級数は普通の意味では収束しない、発散級数であるが、超関数の意味では収束している。また 部分和は普通の関数になっていて、グラフも描ける。
h[n_,x_]:=4/Pi Sum[Cos[k x],{k,1,n,2}]
Manipulate[Plot[h[n,x], {x,-7,7}, PlotRange -> Full, AspectRatio->Automatic, PlotPoints -> 100], {n, 2, 20, 2}]
工学系の本には「デルタ関数列」というものが良く出て来る(数学の解析の本で、軟化作用素を 定義する際の関数列と同じものと考えて良い)。これは n → ∞ の極限で δ となる関数列 {φn} で
ある。 ∫ ∞
−∞
φn(x)dx= 1
という性質を満たしつつ、n が大きくなるにつれ、x= 0 の近くに集中していくものになっている。
例えば工学の本で良く見かけるのは
φn(x) = sinnx πx . このグラフも描いてみて、hn のグラフと見比べてみると良い。
なお、ここまで来ると推察できると思うが、h は g の (超関数の意味での)導関数である。実際、
Heaviside の階段関数
H(x) = {
1 (x≥0) 0 (x <0) の導関数は δ になることが知られているので、
g′(x) = 2δ(x)−2δ(x−π) =h(x) ([−π+ε, π+ε]での等式).
そして、これまでと同様に、h のFourier 級数展開は、g の Fourier級数展開を項別微分したものに なっている。
付 録 J メモ
周期関数の Fourier変換には超関数が必要になるが、遠方で 0 にカットして Fourier 変換したら どうなるだろう？
fm(x) = {
sinx (|x| ≤mπ) 0 (|x|> mπ)
FourierTransform[If[Abs[x] < Pi, Sin[x], 0], x, y]
結果はなるほどである。
f(x) = ei2πxT として、
fm(x) = {
f(x) (|x| ≤mT)
0 (それ以外)
f: R→C が周期 T とする。
cn = 2 T
∫ T /2
−T /2
f(x)e−i2πTxdx とおくと、
f(x) =∑
n∈Z
cnei2πTx.
fm(x) = {
f(x) (|x| ≤mT /2)
0 (それ以外)
とする。
f(ξ) =b 1
√2π
∫ ∞
−∞
f(x)e−iξxdx であるから、
fbm(ξ) = 1
√2π
∫ mT /2
−mT /2
f(x)e−iξxdx 特に
fbm(n) = 1
√2πx
∫ mT /2
−mT /2
f(x)e−inxdx
付 録 K Laplace 変換、 z 変換
Laplace変換、z変換まで含めて話が出来れば良いのだけど…
K.0.1 Laplace 変換
x: [0,∞)→C の Laplace 変換とは
L[x(t)](s) = X(s) =
∫ ∞
0
x(t)e−st dt.
x(t) = 1 2πi
∫ α+i∞ α−i∞
X(s)estds.
ただし α は Laplace 変換の収束域が Res > γ であるとして、α > γ であるように選ぶ (自由度が あるがどれを取っても同じ)。