Modelicaでローレンツ方程式を解く(カオス挙動を確認する)
また前回の更新から大分日が空いてしまいました。 今回はモデリング言語Modelicaを使ってローレンツ方程式を解いていきます。
ローレンツ方程式とは
ローレンツ方程式とは大気変動のメカニズムを簡略して数学モデル化したものです。マサチューセッツ工科大学のローレンツ氏(~2008)によって発表されたのでこの名前がついています。 このローレンツ方程式は下のような式になります。
これをmodelicaによって実装してみます。
modelicaによる実装と結果
ここでは定数をとします。 また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の時間変化は下図のようになります。
シミュレーション開始当初は変動するものの、やがて一定値に収束します。
また、x1に対してx2をプロットした図は下のようになります。
x1,x2ともに一定値へと収束していく様子がわかります。
カオスへの遷移
次に定数rの値とx1,x2,x3の初期値を変えていきます。 次のグラフはr=18として(x1,x2,x3)=(6.7,6.7,17)としたものです。
収束のスピードが前と比べて遅くなっています。
そして次のグラフはr=28として(x1,x2,x3)=(9,8,27)としたものです。 このグラフでは最早(x1,x2)は一定値に収束することはありません。 この独特な形状は「ローレンツ・アトラクター」または「ローレンツ・バタフライ」と呼ばれています。 このように定数と初期値を変更することによって数式モデルの挙動が大幅に変わることが確認できました。
またこのモデルはカオス挙動に特有の初期値鋭敏性を示します。 次のグラフは初期値(x1,x2,x3)=(9.00,8,27)の場合と(9.01,8,27)の場合でのx1の値を比較したものです。
このようにシミュレーション開始時はどちらもほぼ同じ値を示すのですが、途中から2つの挙動が全く違うものとなります。 9.00と9.01というほんの僅かな初期値の違いでこれほどの挙動の違いが生まれるのですから、この系は全く予測できないカオス挙動を示していることになります。
まとめ
Modelicaを使って簡単にローレンツ方程式を解くことができ、カオス挙動を確認することができました。