システムとモデリング

SysML, Matlab/Simulink,他モデリング全般について。たまにプロジェクトマネジメント。

Modelicaでローレンツ方程式を解く(カオス挙動を確認する)

また前回の更新から大分日が空いてしまいました。 今回はモデリング言語Modelicaを使ってローレンツ方程式を解いていきます。

ローレンツ方程式とは

ローレンツ方程式 - Wikipedia

ローレンツ方程式とは大気変動のメカニズムを簡略して数学モデル化したものです。マサチューセッツ工科大学ローレンツ氏(~2008)によって発表されたのでこの名前がついています。 このローレンツ方程式は下のような式になります。

\frac{dx_{1}}{dt}=-\sigma x_{1} +\sigma x_{2}

\frac{dx_{2}}{dt}=-x_{2} +r x_{1}-x_{1}*x_{3}

\frac{dx_{3}}{dt}=-b x_{3} +x_{1}*x_{2}

これをmodelicaによって実装してみます。

modelicaによる実装と結果

ここでは定数を\sigma = 10 , b=8/3  とします。 またr=8, (x1,x2,x3)の初期値は(1,1,1)とします。

model Lorentz
parameter Real sigma=10;
parameter Real b=8/3;
parameter Real r=8"空気層の上部と下部の温度差";
Real x1"対流ロールの回転率";
Real x2"上昇気流と下降気流の温度差";
Real x3"垂直方向の温度分布の線形化からの偏差";

initial equation
x1=1;
x2=1;
x3=1;

equation
der(x1)=-sigma*x1+sigma*x2;
der(x2)=-x2+r*x1-x1*x3;
der(x3)=-b*x3+x1*x2;

end Lorentz;

このときx1,x2,x3の時間変化は下図のようになります。 f:id:Otepipi:20180317115113p:plain

シミュレーション開始当初は変動するものの、やがて一定値に収束します。

また、x1に対してx2をプロットした図は下のようになります。 f:id:Otepipi:20180317115125p:plain

x1,x2ともに一定値へと収束していく様子がわかります。

カオスへの遷移

次に定数rの値とx1,x2,x3の初期値を変えていきます。 次のグラフはr=18として(x1,x2,x3)=(6.7,6.7,17)としたものです。 f:id:Otepipi:20180317115855p:plain

収束のスピードが前と比べて遅くなっています。

そして次のグラフはr=28として(x1,x2,x3)=(9,8,27)としたものです。    f:id:Otepipi:20180317120258p:plain このグラフでは最早(x1,x2)は一定値に収束することはありません。 この独特な形状は「ローレンツ・アトラクター」または「ローレンツ・バタフライ」と呼ばれています。 このように定数と初期値を変更することによって数式モデルの挙動が大幅に変わることが確認できました。

またこのモデルはカオス挙動に特有の初期値鋭敏性を示します。 次のグラフは初期値(x1,x2,x3)=(9.00,8,27)の場合と(9.01,8,27)の場合でのx1の値を比較したものです。 f:id:Otepipi:20180317120741p:plain

このようにシミュレーション開始時はどちらもほぼ同じ値を示すのですが、途中から2つの挙動が全く違うものとなります。 9.00と9.01というほんの僅かな初期値の違いでこれほどの挙動の違いが生まれるのですから、この系は全く予測できないカオス挙動を示していることになります。

まとめ

Modelicaを使って簡単にローレンツ方程式を解くことができ、カオス挙動を確認することができました。