1はじめに

理工学の分野では,連続的に時間変化する自然現象や社会現象等を体系的に理解するために,微分方程式によるモデル化がしばしば用いられる.本研究では,次のような常微分方程式の周期軌道の安定性解析について考える.ここで,x∈RN,f:R×RN→RNはC2級写像であるとする。fがxに関して非線形である場合,方程式(1)の解xの挙動は,パラメータの変化で大きく変化する,すなわち分岐現象を示すことがある[2].この分岐を理解するために,平衡点や周期軌道等,特徴的な解の構造のパラメータに対する変化を調べる.

(1)はベクトル場fが時間に依存する非自律系の方程式である.x(t0)=x0を初期条件とする(1)の解をx(t)=ψ(t-t0)(t0)(x0)と表すとする.周期Tの周期軌道はx(t0+T)=X0を満たすので,周期軌道の条件は次のように表すことができる.周期軌道x(t)の摂動v(t)は,v(t0)=v0であるとき,次のような変分方程式の行列解X0(τ)∈R(N×N)(τ=t-t0)を用いてv(t)=X0(τ)v0と表される.ここで,初期値X0(0)としては通常単位行列Iが用いられる.周期軌道の安定性は,X0(T)の固有値λj(j=1,2,…,N),すなわちFloquet乗数を用いて調べられる.しかし,既存の手法ではFloquet乗数の誤差が大きい場合があることが報告されている[1,3].この問題は陰的Runge-Kutta法等を用いて微分方程式の計算精度を上げても解決できない.

本研究では変分方程式の反復的解法の構成と,変分方程式の初期値をX0(T)の固有ベクトルを用いて設定することで高精度にFloquet乗数を求める手法を提案し,その有効性を数値実験の例によって示す.

2 周期軌道の数値的安定性解析とその問題点

周期軌道の計算方法はすでにいくつか提案されている[2].本研究の手法はshooting法に基づいているので,ここでは比較的簡便かつ広く使われているsimple shooting法について考える.非自律系の場合,周期Tは既知であるので,未知変数はx0である.simple shooting法は,周期解の条件(2)をみたすx0をNewton法等の反復法によって求める.

例として,次式で定義されるMathieu方程式のp=2,20,40における周期軌道を図1(a),(b)に示す.ここで,x=(x,y)T.図1(a)より,周期軌道はpの増加とともに,状態空間の原点に存在する平衡点x*に近づくことが分かる.また,x*付近を通る周期軌道上の解は,x*付近ではゆっくり変化し,それ以外では比較的速く変化している.

次に,Floquet乗数の計算方法を説明する.既存の手法は,(2)を用いて求めた周期軌道を変分方程式(3)に代入し,その変分方程式(3)を初期値問題としてRunge-Kutta法等を用いて数値的に解くことによりX0(T)を求める.さらにX0(T)の固有値を標準的な手法で解くことによりFloquet乗数を求める.

図2(a)に既存の手法で計算したMathieu方程式(4)のFloquet乗数λ1,λ2を示す.Floquet乗数の比|λ1|/|λ2|がpとともに急速に増加し,p>20でλ2が正しく求められていないことが分かる.この例が示すように,既存の手法を用いるとFloquet乗数の誤差が大きくなる場合があることは,他の研究でも報告されている[1,3].

3提案手法

既存のFloquet乗数の計算方法の問題点として,以下の2点が考えられる.

(i)変分方程式(3)の解X0(T)を計算誤差の制御が難しい初期値問題として解いている点

(ii)絶対値が大きな成分を含むX0(T)の固有値を計算する際の行列の変形における丸め誤差

本研究では,これら二つの問題点に対して,それぞれ次のような対策を検討した.

(I)変分方程式の反復的解法の構成

(II)変分方程式の初期条件としてX0(T)の固有ベクトルを用いることによる固有値計算における丸め誤差の低減

図2(b)は,以下で説明する本研究の提案手法で求めたMathieu方程式(4)のFloquet乗数の計算結果を表している.図2(a)と比較すると,既存の手法より提案手法の方が精度の良い結果を与えていることがわかる.

3.1変分方程式の反復的解法

本研究では,常微分方程式の解x(t)=ψ(t-to)(t0)(X0)の次のような性質に着目した.これは力学系の定義としても用いられる一般的な性質である.(5)より次式が得られる.周期軌道の場合はx(t0+T)=X0であるので,x0に関する次のような条件が得られる.変分方程式(3)の解X0(τ)=Ψτ(t0),(x0)(X0(0))も(6)と同じ性質〓を持つため,X0(T)に関する条件が次のように得られる.ここで,0は零行列を表し,X0(0)=Iとする.sの値は自由に選べるが,本研究では簡単のため=T/2とした.条件(7),(8)に対する反復法sを用いることで,X0だけでなくX0(T)の精度も収束条件を用いて制御できる.

3.2 X0(T)の固有ベクトルを用いた変分方程式の初期値の設定

3.1の反復法を用いて精度の良いX0(T)を得ることができるが,X0(T)の成分の絶対値が大きくなると,固有値を計算する段階で生じる丸め誤差のためにFloquet乗数が正しく求まらない場合があった.そこで,以下のような変分方程式の初期値としてXo(T)の固有ベクトルを用いる方法を考えた.

Floque七乗数λj(j=1,2,…,N)に付随するX0(T)の正規固有ベクトルをv(j,0)とすると,〓のように表すことができる.上式の左からv(j,0)の共役ベクトルをかけると次式を得る.ここで,ψT(t0)(x0)(vj,0):=ΨT(t0),(x0)(I)v(j,0), また,〈・,・〉は内積を表す.

(10)は行列の変形を必要としないため,通常の固有値計算で生じる丸め誤差は防げる.したがって,v(j,0)の精度が良ければFloquet乗数λjも高精度に計算できる.そこで,v(j,0)を精度よく求めるために3.1の反復法を以下のように修正した.変分方程式の初期値X0(0)は正則であればいいので,互いに独立なN個の固有ベクトルv(j,0)を用いてX0(0)-v1 ,0…vN,0)を構成する.このとき,X0(T)はX0(0)を用いて次のように表される.ここで,対角行列D0の各(j,j)成分λjは,(10)よりX0とX0(0)の第j列V(j,0)に依存する.したがって,(8)は次のように書き換えられる.以上より,条件(7),(12)を満たす未知変数X0とX0(0),すなわちX0(T)の固有ベクトルv(j,0)を反復法で求め,(10)を用いてFloquet乗数を計算することができる.本研究では,まず条件(7)を満たす周期軌道としてx0を求め,つぎにx0を(12)に代入して/X0(0)を求めた.また,反復法で用いられるshooting法の精度を上げるために,t0

3.3Fioquet乗数の計算精度の指標

Floque七乗数はどの時刻tk(k=0,1,…,M)から一周期を考えても理論的には等しい.しかし,一般に,Xk(T)の各固有値λ(j,k)(j=1,2,…,N)の実際の計算結果は互いに異なる.そこで,本研究ではFloquet乗数λjの精度の指標として,次式R(λj)を用いた。

4数値実験の例

Mathieu方程式(4)は,非自律系であるが,状態空間(x,y)の原点に平衡点を持つ.pの増加とともに平衡点に近づく周期軌道とそのFloque乗数を求めた.図2(a),(b)に従来の手法と提案手法によるFloquet乗数(|λ||≧|λ2|)の計算結果の比較を示す.軌道の分割数はM+1=50とした.提案手法ではp>20でも正確なλ2が得られた.

図3(a),(b)はFloquet乗数の精度の指標R(λj)の比較を示す.従来の手法ではpの増加とともにR(λ2)が急速に増加したが,提案手法ではR(λ1)とR(λ2)のどちらもほぼ0に抑えられた.これは,提案手法の有効性を示している.

また,提案手法はわずかな修正によって,ベクトル場fが時間に依存しない自律系に対しても適用することができる.本論文では,3次元自律系のFitzHugh/Nagumo方程式に対しても数値実験を行い,提案手法によって高精度なFloquet乗数が得られることを確認した.

提案手法の計算時間は,どちらの計算例に対しても既存の手法の1.5倍未満だった.計算量に対する詳細な考察は行っていないが,今回行った数値実験の範囲では,実用的な計算時間で高精度なFloquet乗数が得られた.

5まとめ

本研究では,常微分方程式の周期軌道の数値的安定性解析について考えた.周期軌道の安定性は,変分方程式の解を用いて定義されるFloquet乗数を用いて調べることができる.Floquet乗数の計算方法はすでにいくつか提案されているが,本研究ではそれらの問題点を整理し,改良手法を提案した.提案手法のアイディアは,変分方程式に対する反復法の構成とFloquet乗数に付随する固有ベクトルを用いた固有値計算である.数値実験により,既存の方法では正確に求めることができなかったFloquet定数を,提案手法により実用的な計算時間で精度よく求めることができることがわかった.

図1 提案手法によるMathieu方程式(4)の周期軌道の計算結果.T:周期.p:方程式のパラメータ.

図2 Mathieu方程式(4)のFloquet乗数λj(j=1,2)の計算結果の比較.p:方程式のパラメータ.

図3 Mathieu方程式(4)のFloquet乗数λj(j=1,2)の計算精度の比較.p:方程式のパラメータ.R(λj):(15).

本論文は、「Accurate Numerical Computation for Floquet Multipliers of Periodic Orbits of Ordinary Differential Equations常微分方程式の周期軌道とFloquet乗数の高精度数値計算法と題し、6章と付録から構成されている。理工学分野において、連続的に時間変化する自然現象や社会現象などを解析するために、常微分方程式によるモデル化がしばしば用いられている。しかし、非線形常微分方程式の場合、パラメータの変化により解の挙動は大きく変化し、解の分岐現象が発生する場合がある。複雑な分岐現象の構造を理解するためには、基本的な定常状態である平衡点と周期軌道の安定性を調べる必要がある。特に周期軌道の場合は、安定性を決定するFloquet乗数を求めなければならないが、その数値計算は困難な場合があることが知られている。本論文は、従来の手法では誤差の大きい計算結果しか得られない場合に対しても、Floquet乗数を精度よくかつ効率よく求めることができる新しい数値計算法を提案し、その有用性を明らかにしている。

第1章「Introduction」では、力学系における安定性解析に対する研究の背景と目的を述べ、従来研究に対する本論文の位置付けを与えている。また、本論文の構成を示している。

第2章「Numerical stability analysis of periodic orbits of ordinary differential equations」では、Floquet乗数の既存の計算手法とその問題 点について述べている。既存の手法では、変分方程式の初期値問題を数値的に解いて行列解を求め、その解の固有値としてFloquet乗数を計算している。しかし、平衡点付近を通る周期軌道に対しては行列解の計算精度の低下と固有値計算の丸め誤差が原因となり、Floquet乗数の計算精度が悪くなる欠点があることを明らかにしている。

第3章「A new iterative method for periodic orbits and Floquet multiplers」では、次の2つのアイデアに基づくFloquet乗数の新たな計算手法を提案している。すなわち、(1)変分方程式の初期値問題として解くことによる誤差の蓄積を回避するために、変分方程式の解に対する新たな反復計算法を構成することと、(2)行列解の成分の絶対値が大きい場合の固有値計算の丸め誤差を回避するために、行列解の固有ベクトルを変分方程式の初期値に用いることである。

第4章「Numerical examples」では、非線形Mathieu方程式とFitzHugh-Nagumo方程式を具体的に取り上げ、それらの方程式に対する数値計算結果を示すことにより、既存の計算手法に比べて、第3章で提案した計算手法が、実用的な計算時間で高精度にFloquet乗数を計算できることを実証している。

第5章「Discussions」では、第3章で提案した2つのアイデアの重要性について考察し、2つのアイデアのそれぞれが、Floquet乗数を高精度に求めるために重要であることを、計算例を用いて明らかにしている。

第6章「Conclusions」では、本論文の成果をまとめると共に、今後の研究課題についても考察をしている。

付録(Appendix)では、論文の自己完結性をよくするために、周期的Schur分解の手法を要約し、紹介している。

なお、本論文の成果は、村重淳との共同研究であるが、論文提出者が主体となって新しい数値計算手法の提案、解析、およびシミュレーションを行なったものであり、論文提出者の寄与が十分であると判断する。

したがって、博士(科学)の学位を授与できると認める。