第 5 章 周期解の近似計算 86
6.2 占部・ガレルキン法と調和平衡法
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0
0.1 0.2 0.3 0.4 0.5 0.6 0.7
0
1
3
C
b
ζ
図6.4 平衡点の分岐図．曲線上で平衡点対の発生・消滅が起こる．
とし，各引数について連続でかつ必要な回数だけ微分可能な性質を持つと仮定しよう．また，記法を簡 単にするため方程式
F(t, x,x) =˙ f(t, x)−x˙ (6.21)
を考える．
さて，式(6.21) が周期2π の周期解x(t) を持っていたとすると，これは次式のようにフーリエ級数
に展開できる．
x(t) =ϕ(t) =a0+ X∞
k=1
(akcoskt+bksinkt) (6.22) 占部・ガレルキン法とは，式 (6.22)の係数ベクトル (a, b) = (a0, a1, b1, . . .) を有限個近似計算する手 法である．
式 (6.22)を式 (6.21)の左辺に代入し，フーリエ級数展開した式を
F(t, ϕ(t),ϕ(t)) =˙ Fc0+ X∞
k=1
(Fckcoskt+Fsksinkt) (6.23)
とすれば，式(6.21) より三角関数の各係数は零となる．
Fc0(a, b) = 1 2π
Z 2π
0
F(τ, ϕ(τ),ϕ(τ˙ ))dτ = 0 Fck(a, b) = 1
π Z 2π
0
F(τ, ϕ(τ),ϕ(τ˙ )) coskτ dτ = 0 k= 1,2, . . . Fsk(a, b) = 1
π Z 2π
0
F(τ, ϕ(τ),ϕ(τ˙ )) sinkτ dτ = 0
(6.24)
そこで，この可算無限個の方程式を有限個で近似し，式 (6.24)の近似解を得る方法が，占部によって研 究されたガレルキン法である．また彼は，逆に得られた近似解から式 (6.24) の真の解の存在の言える ことも示した．以下計算の手続きのみを示しておこう．ここでもニュートン法を用いる．その際の計算 がすべて数値的になされることを見ておきたい．
周期解 (6.22)を m 次の高調波で打ち切った近似周期解を
x(t) =ϕ(t) =a0+ Xm
k=1
(akcoskt+bksinkt) (6.25) とする．これを m 次のガレルキン近似解という．これを式(6.21)に代入し， m 次の高調波までの係 数を零とおくと，式 (6.24)と同様な(2m+ 1)n 個の方程式
Fc0(a, b) = 1 2π
Z 2π
0
F(τ, ϕ(τ),ϕ(τ˙ ))dτ =Pc0(F) = 0 Fck(a, b) = 1
π Z 2π
0
F(τ, ϕ(τ),ϕ(τ˙ )) coskτ dτ =Pck(F) = 0 k= 1,2, . . . , m Fsk(a, b) = 1
π Z 2π
0
F(τ, ϕ(τ),ϕ(τ˙ )) sinkτ dτ =Psk(F) = 0
(6.26)
を得る．これをm 次の決定方程式と呼ぶ．ただし
(a, b) = (a0, a1, b1, . . . , am, bm) (6.27) となっている．また，各係数を取り出す積分の操作を Pc0, Pck, Psk, . . . と略記した．
そこで，式(6.27)の係数(a, b)を未知変数と考えて式(6.26)を解く．ニュートン法で解くことにし，
ヤコビ行列の各要素がどうなるか見てみよう．まず
∂F
∂a0 = ∂Fc0
∂a0 + Xm
k=1
µ∂Fck
∂a0 coskt+ ∂Fsk
∂a0 sinkt
¶
= ∂F
∂x
∂x
∂a0 +∂F
∂x˙
∂x˙
∂a0 = ∂F
∂x
∂F
∂aj = ∂Fc0
∂aj + Xm
k=1
µ∂Fck
∂aj coskt+ ∂Fsk
∂aj sinkt
¶
= ∂F
∂x cosjt−∂F
∂x˙jsinjt
∂F
∂bj
= ∂Fc0
∂bj
+ Xm
k=1
µ∂Fck
∂bj
coskt+ ∂Fsk
∂bj
sinkt
¶
= ∂F
∂x sinjt+ ∂F
∂x˙jcosjt
(6.28)
の関係のあることに注目して，各要素は次式となる．
∂Fc0
∂a0 = Pc0
µ∂F
∂x
¶
∂Fck
∂a0 = Pck
µ∂F
∂x
¶
, ∂Fsk
∂a0 =Psk
µ∂F
∂x
¶
, k= 1,2, . . . , m
∂Fc0
∂aj = Pc0
µ∂F
∂x cosjt−∂F
∂x˙jsinjt
¶
, j= 1,2, . . . , m
∂Fc0
∂bj = Pc0
µ∂F
∂x sinjt+∂F
∂x˙ jcosjt
¶ ,
∂Fck
∂aj
= Pck µ∂F
∂x cosjt− ∂F
∂x˙jsinjt
¶
, j, k= 1,2, . . . , m
∂Fsk
∂aj
= Psk µ∂F
∂x cosjt−∂F
∂x˙jsinjt
¶ ,
∂Fck
∂bj
= Pck µ∂F
∂x sinjt+∂F
∂x˙jcosjt
¶ ,
∂Fsk
∂bj = Psk µ∂F
∂x sinjt+∂F
∂x˙ jcosjt
¶
(6.29)
したがって，これらは近似解(6.25)の係数 (6.27) を与えると F(t, ϕm(t),ϕ˙m(t)), ∂F
∂x(t, ϕm(t),ϕ˙m(t)), ∂F
∂x˙(t, ϕm(t),ϕ˙m(t)) の計算とフーリエ係数を求めることで遂行できる．
【例 6.2】占部・ガレルキン法による周期解の計算例 ダフィング方程式
¨
x+ 0.1 ˙x+x3= 0.3 cost (6.30)
の周期2π の周期解を計算してみよう．ここでは計算の手順を理解するためにガレルキンの低次の近似 解を実際に代入しながら決定方程式を導くことにする．
いま，x(t)を式(6.30)の周期解とすると，この式の非線形性が奇関数であることから，−x(t+π)も
周期解である．特に x(t) =−x(t+π) の性質を持つ周期解が存在する．この性質を持つ周期解はフー リエ級数に展開すると次のように奇数次の高調波しか含まれない．
x(t) = X∞
k=0
©a2k+1cos(2k+ 1)t+b2k+1sin(2k+ 1)tª
(6.31) そこで，すこし粗い近似解であるが 2次のガレルキン近似解：
x(t) =a1cost+b1sint+a3cos 3t+b3sin 3t (6.32) を求めてみよう．式 (6.32) を式 (6.30) に代入して cost,sint,cos 3t,sin 3t の係数のみを両辺比較し て，これらの係数を零とおくと2 次の決定方程式を得る．
Fc1(a1, b1, a3, b3) = −a1+ 0.1b1+ 3 4 h
a1(r12+ 2r23) +a3(a21−b21) + 2b3a1b1 i
−0.3 = 0 Fs1(a1, b1, a3, b3) = −b1−0.1a1+ 3
4 h
b1(r21+ 2r23) +b3(a21−b21)−2a3a1b1 i
= 0 Fc3(a1, b1, a3, b3) = −9a3+ 0.3b3+1
4 h
a1(a21−3b21) + 3a3(2r12+r32) i
= 0 Fs3(a1, b1, a3, b3) = −9b3−0.3a3+1
4 h
b1(3a21−b21) + 3b3(2r21+r23) i
= 0
(6.33) ここに，
r21=a21+b21, r23=a23+b23 とおいた．式(6.31) をニュートン法で解くと次の3個の解を得る．
解1 : (a1, b1, a3, b3) = (−0.9023, 0.3047, −0.0149, 0.0240) 解2 : (a1, b1, a3, b3) = (−0.3219, 0.0349, −0.0009, 0.0003) 解3 : (a1, b1, a3, b3) = (1.1198, 0.5246, 0.0209, 0.0673)
これらの周期解を例 4.7（計算法については例 6.6 参照）と比較するとかなりよい近似解が得られてい ることが分かる．実際の占部・ガレルキン法ではすべての手順が数値計算で処理できるので決定方程式 やそのヤコビ行列を陽に導出する必要はない．したがって適当な m 次の近似解を求めて，次に次数m を大きくしていって最高次数の係数が十分に小さくなった時点で近似の次数を定めるとよいであろう．
■
6.2.2 調和平衡法と形式的平均化法
占部・ガレルキン近似解を最初の基本周波数のみ，あるいは顕著な周波数の成分のみの近似解で打 ち切った場合は，すべてが簡単となる．決定方程式も直接導出できることが多い．これを調和平衡法 (harmonic balance method) あるいは等価線形化法(equivalent linearization) という．このことを具 体例でみてみよう．
【例 6.3】調和平衡法によるダフィング方程式の周期解 ダフィング方程式
¨
x+kx˙+c1x+c3x3=Bcost (6.34) を考えよう．いまこの方程式の周期解として
x(t) =ucost+vsint (6.35)
を仮定する．これを式 (6.34)に代入してcostとsint の係数をそれぞれ零とおくと次式を得る．
Fc(u, v) = µ
−1 +c1+3 4c3r2
¶
u−kv−B = 0 Fs(u, v) = ku+
µ
−1 +c1+3 4c3r2
¶ v= 0
(6.36)
ここに，
r2=u2+v2
とおいた．これは摂動法で得た振幅特性 (5.59) と類似の式となっている．したがって解析は，例 6.1 と同様に行えばよい． ■
【例 6.4】形式的な平均化法
式 (6.34)について，求めた周期解 (6.35)の近傍で解の様子をみるため，今度は
x(t) =u(t) cost+v(t) sint (6.37) と考えて，係数が変化するとしよう．式 (6.37)の時間微分を
˙
x = u(t) cos˙ t+ ˙v(t) sint−u(t) sint+v(t) cost
¨
x = −2 ˙u(t) cost+ 2 ˙v(t) sint−u(t) cost−v(t) sint
と考え，これを式 (6.34)に代入する．ここでu(t), v(t) はゆっくり変化すると考えたので，
¨
u, v , k¨ u, k˙ v,˙ sin 3t, cos 3t
などに関係した項はすべて無視した．そこでcostとsint の係数を比較して
˙
u = 1
2 (
−ku− µ
1−c1− 3 4c3r2
¶ v
)
˙
v = 1
2 ( µ
1−c1−3 4c3r2
¶
u−kv+B
) (6.38)
を得る．ただし
r2=u2+v2
とおいた．ここでも平均化法と同じ結果の得られることに注意しよう．ただ，このような導出過程には 数学的根拠がはっきりしていないため，ここでは「形式的」平均化法と呼ぶことにした． ■
調和平衡法や形式的な平均化法は，非線形系に調和振動に近い，すなわちあまり高調波を含んでいな い振動が見られる場合，簡便で実用的な解析方法として用いられている．数学的に周期解の存在や近似 の正当性を証明できないが，適切な仮定の下で，実際の系の挙動を十分に近似する場合がある．どのよ うな振動の成分を仮定し，解析を進めるかは目前の現象をみて適切に定めなければならない．