1分岐集合の数値計算
吉田勝俊

2007 年 12 月 11 日

2 自律系の平衡点

1 力学系一定の規則に従って時間変動するシステムを，力学系自律系 (1) を考える．そこに至ると時間変動が止まるという．力学系の配位 (物体なら位置と向き) を表すのに相空間の点を，平衡点または不動点という．この定義よ必要な変数の組 (x1 , x2 , · · · , xn ) ∈ Rn を，一般化座標とり，平衡点では一般化座標の時間微分が全て 0 になる．

いう．一般化座標が張る数ベクトル空間を，配位空間と

ẋ1 = ẋ2 = · · · = ẋn = 0いう．一般化座標の時間微分を加味した，これを自律系 (1) に代入すれば，平衡点が満すべき条件(x1 , ẋ1 , x2 , ẋ2 · · · , xn , ẋn ) ∈ R2nとして，

f1 (x1 , x2 , · · · , xn ) = 0f2 (x1 , x2 , · · · , xn ) = 0..
.

fn (x1 , x2 , · · · , xn ) = 0が張る数ベクトル空間を，相空間という．
時間変動が，n 連立の微分方程式，

ẋ1 = f1 (x1 , x2 , · · · , xn )ẋ2 = f2 (x1 , x2 , · · · , xn )..
.
ẋn = fn (x1 , x2 , · · · , xn )(1)が得られる．
次の 2 次元力学系を例題に，数値計算法を示そう．に支配される力学系を考える．解は相空間上の曲線とな{
ẋ1 = x2 ,ẋ2 = −x2 − px1 − x31 + F.

るが，これを相軌道という．ある瞬間の力学系の状態はこの系を，3 次の Duﬃng 系という．一般に，復元力が多相軌道上の 1 点に対応する．このような，右辺が時間 t項式で表わされる振動系を Duﬃng 系という．

(3)を陽には含まない力学系を，自律系という．
これに対して，右辺が時刻 t を陽に含む系：

ẋ1 = f1 (x1 , x2 , · · · , xn , t)ẋ2 = f2 (x1 , x2 , · · · , xn , t)..
.
ẋn = fn (x1 , x2 , · · · , xn , t)

2.1 平衡点の探索
自律系 (3) の平衡点は，停止条件 ẋ1 = ẋ2 = 0 より，(2)x2 = 0,−px1 − x31 + F = 0を満足する．簡単のため p = −1, F = 0 としよう．
を，非自律系という．具体例を挙げると，ẋ1 = x2 ,ẋ2 = −cx2 − kx1 − x31は自律系である．その一方で，

ẋ1 = x2 ,このとき，解くべき非線形方程式は，

ẋ2 = −cx2 − kx1 − ϵx31 + A cos ωtx1 − x31 = 0となるから，この場合の平衡点は手計算できる．

x1 (1 − x21 ) = 0∴ x1 = 0, ±1

以上，この系には平衡点が 3 個ある．

とか，ẋ1 = x2 ,ẋ2 = −cx2 − kx1 cos ωt

この例題は幸運にも手で解けたが，一般の非線形方程式は手では解けないので，数値計算の実例を示していく．などは陽に t を含むので非自律系である．ようするに，■数値積分の欠点

非自律系には，xi や ẋi と全く無関係な時間変動が，外部あるから，元の微分方程式 (3) の数値解を求めて観察すから注入される．

れば，いずれ平衡点に収束しそうである．

平衡点とは，現象的には定常状態でDRAFT@2007/12/11

function dx = eqn(x, t)
global P;
dx(1) = x(2);dx(2) = -x(2) + x(1) - x(1)**3;
endfunction
t=linspace( 0.0, 20.0, 200 );initp=linspace( -2.0, 2.0, 10 );all=[];
for i=1:10
sol = lsode( "eqn", [initp(i); 0.0], t );all=[all, sol(:,1)];
endforxlabel("Time");ylabel("x1");
plot(t,all);x1 軸との交点は，やはり 3 点ある．数値で求めるにはF (x) が十分小さくなる x を拾ってやればよかろう．
function y = F(x)y = -(x - x.**3);
endfunction
n=5000; eps=1e-3;x=linspace( -2.0, 2.0, n );fp0=[];
for i=1:n
if( abs(F(x(i))) < eps )printf( "%e %e\n", x(i), F(x(i)) );fp0 = [fp0, x(i)];endifendforsave fp0.dat fp0

このように，この方法の欠点として，3 つある平衡点のこれを実行すると，(必要なら n や eps を調整しながら)うち x = ±1 は検出できるが，x = 0 が検出できない．理-9.998000e-01 3.999600e-04-4.000800e-04 4.000800e-04
4.000800e-04 -4.000800e-04
9.998000e-01 -3.999600e-04

由は後述するように，x = ±1 が安定平衡点，x = 0 が不安定平衡点だからである．不安定な平衡点 x = 0 は，自然状態では発現しない*1 ．自由状態の振り子を例にとると，下死点が安定平衡点，上死点が不安定平衡点である．のような出力を得る．計算結果を収めたベクトル変数

数値積分法は，このような物理現象を忠実にシュミレー

fp0 の内容は fp0.dat というファイルに保存される．

トしてしまうので，不安定平衡点の数値解もまた，通常x = −1, 0, 1 付近の点が出力されており，粗いながらも，の観測には掛からない*2■作図による方法

平衡点の推定値が得られている．

そこで少々泥臭いが，平衡点の方程作図による推定値の精度は |F (x)| <10−4 程度だが，これらを「磨く」ために，ニュートン法式のグラフを書く方法を紹介しておく．ようするに，(速■ニュートン法度の項を零においた) 静的な復元力，を使おう．Octave にはそれ用に fsolve という関数が用

F (x1 ) = −(x1 − x31 )意されている．のグラフを書いて*3 ，x1 軸との交点を求めれば，それがload fp0.dat
dn=length(fp0);
function y = F(x)y = -(x - x.**3);
endfunction
fp=[];
for i=1:dn
fp(i)=fsolve("F", fp0(i));
printf( "%e %e\n", fp(i), F(fp(i)) );endforsave fp.dat fp

平衡点である．
function y = F(x)y = -(x - x.**3);
endfunction
x=linspace( -2.0, 2.0, 200 );xlabel("x"); ylabel("F(x)"); grid;plot(x,F(x));*1 ただし，クーロン摩擦などがあるときは，不安定平衡点に静止す

これを実行すると，る状態が安定に観察できる．
*2 不安定平衡点が求まる数値計算法も存在するが (Shooting 法な-1.000000e+00 -0.000000e+00
0.000000e+00 -0.000000e+00
0.000000e+00 -0.000000e+00
1.000000e+00 -0.000000e+00

ど)，プログラムの微調整に熟練が必要である．
*3 これを復元力とみたときに辻褄が合うように，マイナスを付けたが，x1 軸との交点に影響はない．

2DRAFT@2007/12/11となり，計算機精度の範囲で正確な平衡点が求まる．ことなるが，下線部は元の式 (1) だからキャンセルして，ズ

こでは幸運にも，厳密解と同じ数値が得られたが，通常レだけの微分方程式，は 10−9 程度の精度なら満足することになる．
■まとめ∂fi∂fi
ξ1 + · · · +ξnξ˙i =∂x1∂xn

最初からニュートン法を使えばよさそうなものだが，そうできない理由がある．ニュートン法は，平が得られる．i = 1, · · · , n について全部並べて，行列表衡点の精度を「磨く」用途には適しているが，平衡点を記すると，

 ˙   ∂f1
ξ1
∂x1
∂f2
 ξ2  
∂x1
  

 ..  =  .
 .   ..

見付ける用途には適していないのである．これはよく知られたニュートン法の欠点で，ニュートン法の初期値，fp(i)=fsolve("F", fp0(i));^^^^^^ 初期値

ξn|を，厳密解から離れた値にセットすると，ニュートン法∂fn∂x1∂fn∂x2{z ∂f1
ξ1∂xn∂f2   
∂xn   ξ2 .
..  .   .. ∂fn
ξn∂xn}

行列部分 A を，自律系 (1) のヤコビ行列という．

理想だが，現状では，数値積分や作図，またはその他の自律系 (1) を，太字のベクトル表記，アイデアで，臨機応変に平衡点を見付けていくしかない．

 x1
 x :=  ...  ,特に数値積分が不安定平衡点を見逃すことには，十分留意すべきである．いったん見付かりさえすれば少々精度が悪くても，ニュートン法で磨ける．


f1 (x1 , · · · , xn )..
f (x) := .
fn (x1 , · · · , xn )x2によって，ẋ = f (x) と書くとき，ヤコビ行列 A は2.2 平衡点の安定性

とか Dx f と表記されることが多い*4 ．

y = f (x1 , · · · , xn ) という関数関係におい∂f∂x■線形化系の安定性

試しに Duﬃng 系，{ẋ1 = x2 ,ẋ2 = −x2 + x1 − x31 .て，xi を独立変数，y を従属変数と呼んだ．ここで，独立変数 xi を一斉に，ちょっとだけずらそう．

(i = 1, 2, · · · , n)(7)となるが，これを，自律系 (1) の第 1 変分方程式という．

全ての平衡点を漏れなく探索できる数値解法があれば

xi 7→ xi + ξi..
.

···
···
..
.
···

A

は厳密解に到達できない．

■第 1 変分

∂f1
∂x2
∂f2
∂x2(4)

ξi ≪ 1 は微小な変数である．このとき，特別な場合を除の第 1 変分方程式を求めてみよう．P は定数とする．式

けば，従属変数 y もズレるはずだ．

(7) に当てはめると，ヤコビ行列は，[] [従属変数 y に起こるズレ η を，ごく自然に，

η := f (x1 + ξ1 , · · · , xn + ξn ) − f (x1 , · · · , xn )∂(第 1 式の右辺)∂x1∂(第 2 式の右辺)∂x1(5)∂(第 1 式の右辺)∂x2∂(第 2 式の右辺)∂x2=

0
1 − 3x2]
1(8)−1となり，第 1 変分方程式は，と定義する．ここで，独立変数のズレ ξi が十分に微小で，[˙] [ξ10=
ξ2
1 − 3x2関数 f (x1 , · · · , xn ) が全微分可能ならば，従属変数のズレ η は，独立変数のズレ ξi の一次式によって，次のよう

1−1][ ]
ξ1ξ2(9)に近似計算できる．

η≈∂f∂fξ1 + · · · +ξn∂x1∂xnとなる．ようするに，Duﬃng 系に与えるズレ ξ が十分に(6a)小さければ，ズレ ξ の時間発展は線形力学系 (9) で近似1 次近似式 (6a) を，y = f (x1 , · · · , xn ) の第 1 変分といされる．このことから，第 1 変分方程式を線形化方程式もしくは線形化系と呼ぶことがある．

う．定義式 (5) を用いれば，近似式 (6a) を，線形力学系の安定性は，システム行列*5 の固有値で判

f (x1 + ξ1 , · · · , xn + ξn )∂f∂f
ξ1 + · · · +ξn≈ f (x1 , · · · , xn ) +∂x1∂xn

別できた．Duﬃng 系の線形化系 (9) のシステム行列はヤコビ行列 (8) だが，その成分には，元の Duﬃng 系の(6b)変位 x が含まれている．ズレ ξ の安定性は，x に何を代と表記しておくこともできる．
■ヤコビ行列入するかで変化する．
この問題について，次の定理が知られている*6 ．

自律系 (1) にズレ (4) を代入すると，(xi + ξi )˙ = fi (x1 + ξ1 , · · · , xn + ξn )=⇒ ẋi + ξ˙i = pf (x1 , · · · , xn )+*4 ∂f は式 (7) に出てきた行列の略記である．ベクトルでベクトル∂xを偏微分する方法を独自にあみだす必要はない．
*5 ẋ = Ax の A のこと．
*6 線形化系の固有値実部が零のときには，分岐という非線形現象が∂fi∂fiξ1 + · · · +ξn ∵ (6b)∂x1∂xn起こり，定理が成立しなくなるのだが，詳しくは後述する．

3DRAFT@2007/12/11• 非線形力学系と線形化系の解軌道が似るのは，平衡点 x = x̄ の近くである．

2.3 平衡点の接続

この定理を前提に，平衡点まわりの線形化系の固有値を■平衡点の追跡求める．

パラメータ p を少しずつ変えながらニュートン法を繰り返すと，平衡点 x をパラメータの関load fp.dat #fp
fixp=fp; fp=[]; fp=fixp([1,2,4]);function A = jac( x )A = [0, 1; 1-3*x**2, -1];
endfunctionev(1,:)=eig(jac(fp(1)))’;ev(2,:)=eig(jac(fp(2)))’;
ev(3,:)=eig(jac(fp(3)))’;
fp’ev
save fp_ev.dat fp ev

数 x(p) として描ける．得られたグラフを分岐図という．
load fp.dat
global x p; p=-1;function y = F(x,p)
y = p*x + x.**3;
endfunction
function y = F_x(x)global p;
y = F(x,p);
endfunction
parn = 101;
bifpar(1,:)=linspace(-1,-0.1,parn); # px0=fp(1); bif(1,:)=[]; # resultfor i=1:parnp = bifpar(1,i); x0=fsolve("F_x", x0);
bif(1,i) = x0;
endforplot( bifpar(1,:), bif(1,:) );function y = F_p(p)global x;
y = F(x,p);
endfunction
p0 = bifpar(1,parn); # last paramx0 = bif(1,parn);# last fixp
bifpar(2,:)=linspace( x0, x0+1, parn ); # xbif(2,:)=[]; # resultfor i=1:parnx = bifpar(2,i); p0=fsolve("F_p", p0);bif(2,i) = p0;
endforplot( bifpar(1,:), bif(1,:), bif(2,:), bifpar(2,:) );p0 = bif(2,parn);# last paramx0 = bifpar(2,parn); # last fixp
bifpar(3,:)=linspace(p0,-1,parn); # pbif(3,:)=[]; # resultfor i=1:parnp = bifpar(3,i); x0=fsolve("F_x", x0);
bif(3,i) = x0;
endforplot( bifpar(1,:), bif(1,:), bif(2,:), bifpar(2,:),\
bifpar(3,:), bif(3,:) );input("hit any key");
fpset1(:,1) = [bifpar(1,:),bif(2,:),bifpar(3,:)]’;
fpset1(:,2) = [bif(1,:),bifpar(2,:),bif(3,:)]’;
save fpset1.dat fpset1plot( fpset1(:,1), fpset1(:,2) ); pause;

実行すると，> octave -q DufFixpEig.m
ans =-101ev =-0.50000 - 1.32288i -0.50000 + 1.32288i-1.61803 + 0.00000i
0.61803 + 0.00000i-0.50000 - 1.32288i -0.50000 + 1.32288iとなる．平衡点 x = −1, 1 は固有値 −0.50000 ± 1.32288i
より安定フォーカス，平衡点 x = 0 は固有値 −1.61803,0.61803 よりサドルであることが分かる．以下，元のDuﬃng 系の解軌道を示す．
function dx = eqn(x, t)dx(1) = x(2);dx(2) = -x(2) + x(1) - x(1)**3;
endfunction
t=linspace( 0.0, 20.0, 200 );ip1=linspace( -0.6, 0.3, 10 );ip2=linspace( -0.3, 0.6, 10 );x=[]; dx=[];for i=1:10
sol = lsode( "eqn", [ip1(i); 0.4], t );x=[x, sol(:,1)]; dx=[dx, sol(:,2)];sol = lsode( "eqn", [ip2(i); -0.4], t );x=[x, sol(:,1)]; dx=[dx, sol(:,2)];
endforxlabel("x");ylabel("dot x");plot(x,dx); pause;

4

DRAFT@2007/12/11

グラフが垂直に立ってくると追跡困難になるので，途中xlabel("p"); ylabel("x");plot( s1(1,:), s1(2,:), "1",\s2(1,:), s2(2,:), "1",\u2(1,:), u2(2,:), "2" ); pause;※ u1 は空なので，プロットしようとするとエラーになる．

で掃引方向を変更し，平衡点 x を少しずつ変えながらパラメータ p を求めて，関数 p(x) を描いている．

Exercise 1 平衡点 x = 0 について，p = −1 から 1 まで接続せよ．
プログラム例を示す．
load fp.dat
global x p; p=-1;function y = F(x,p)
y = p*x + x.**3;
endfunction
function y = F_x(x)global p;
y = F(x,p);
endfunction
parn = 301;
bifpar(1,:)=linspace(-1,1,parn); # px0=fp(2); bif(1,:)=[]; # resultfor i=1:parnp = bifpar(1,i); x0=fsolve("F_x", x0);
bif(1,i) = x0;
endforplot( bifpar(1,:), bif(1,:) );input("hit any key");
fpset2(:,1) = bifpar(1,:)’;
fpset2(:,2) = bif(1,:)’;
load fpset1.datplot( fpset1(:,1), fpset1(:,2), \fpset2(:,1), fpset2(:,2) );save fpset2.dat fpset2
pause;■平衡点の安定判別■ピッチフォーク分岐 以上の結果を，p = 1 から p =−1 の方向にたどると，次のことがいえる．

• p > 0 では安定平衡点 x = 0 のみが存在する．
• p = 0 において，x = 0 が不安定化し，x = 0 の上下に新たな安定平衡点が生じる．
このような変化の仕方を，超臨界ピッチフォーク分岐(super critical pitchfork bifurcation) という．
一般に，平衡点の性質や個数が変化する現象を分岐現象といい，分岐が起こる境界点を分岐点という．以上の固有値実部が負となる安定平衡点

例では，p = 0 が分岐点である．と，それ以外の不安定平衡点を分離して，色分けしてプロットする方法を示す．

■サドル・ノード分岐

load fpset1.datload fpset2.dat
function A = jac( x, p )A = [0, 1; -p-3*x**2, -1];
endfunction
function [s, u] = stab(fpset)
x_s = x_u = p_s = p_u =[];
for i=1:rows(fpset)
p = fpset(i,1);
x = fpset(i,2);
ev_real = real(eig(jac(x,p)));if ( ev_real(1) < 0 )p_s=[p_s, p]; x_s=[x_s, x];else
p_u=[p_u, p]; x_u=[x_u, x];
endif
endfors=[]; u=[];if ( length(p_u) > 0 )u=[u, [p_u;x_u]];endif
if ( length(p_s) > 0 )s=[s, [p_s;x_s]];
endif
endfunction

これまでの例では，以下の Duﬀ-ing 系において F = 0 としてきた．
{ẋ1 = x2 ,ẋ2 = −x2 − px1 − x31 + F.
Exercise 2 F ̸= 0 のときに分岐図がどう変化するか，確認せよ．

このような F ̸= 0 のときに生じる分岐現象を，サド[s1,u1]=stab(fpset1);[s2,u2]=stab(fpset2);
fpset_s=[s2,s1]; fpset_u=[u1,u2];
grid on;

ル・ノード分岐という．サドルノード分岐点では，サド
ルと安定ノードの組が消失する．

5

DRAFT@2007/12/11

2.4 分岐点の接続■分岐点の判定

各平衡点 x̄ のヤコビ行列 Df (x̄) の固有値 λ:P (λ) := det(Df (x̄) − λI) = 0を計算すると，分岐図が垂直に切り立つ点で実部が零 0になる．これを分岐点という．
分岐点の特殊な場合として，虚部も同時に消えて固有値が 0 になるときは，行列式の性質によって，P (0) = det(Df (x̄)) = 0が成立する．この条件に該当する分岐点として，ピッチフォーク分岐点やサドル・ノード分岐点がある．
load fpset1.dat; #DufCont1.m
load fpset2.dat;
fps = fpset1;function A=jaco( x, p )A=[0, 1; -p-3*x**2, -1];
endfunction
p=fps(:,1); x=fps(:,2);
det_jac=[];for i=1:rows(x)det_jac(i) = det(jaco(x(i),p(i)));
endforylabel("Determinant P(0)");xlabel("x"); plot( x, det_jac ); pause;
xlabel("p"); plot( p, det_jac ); pause;■分岐点の接続

平衡点の方程式と P (0) = 0 を連立して解く．

6