今北産業ローレンツアトラクタ

chouchoが舞えば気象兵器が儲かる
ローレンツ方程式
ほれ3行

'(ダッシュ・プライム)は微分を意味するので、この3行の方程式は連立微分方程式です。時間tで微分します。
なお、非線形ですので行列は帰って、どうぞ。

x,y,zが求められるべき関数で、x(t)、y(t)、z(t)らしいです。

p,r,bはパラメータです。
p=10、r=28、b=8/3っていうのが割りと注目されるらしいです。

今回数値計算でシミュレーションしましたが、初期値はx=10、y=12、z=15を参考にしました。

一度オイラー法での差分方程式によるシミュレーションを試しましたがすぐに変な値に収束してしまうので、ルンゲ・クッタ法でやり直したらExcelでもうまくいきました。陰的とかブッチャー配列とか知りません。古典的でやりました。



ルンゲクッタ法は以上のような手続きを取るそうです。
積分における短冊・台形近似をシンプソンの公式にしたやつの微分方程式版と言われて言葉ではなんとなく把握したけども、いまいち理解していません。

ですので、応用が利かず、たとえばy'=-xz+rx-yの場合任意の関数fはf=-xz+rx-yなのでf(x,z,y,t)なんでしょうが今はx,zのことは考えなくていいとして
yにxとかzとか掛け算されてたらどうしよう！？ｇｋｂｒとか思ってたんですが3本の式とも杞憂に終わってしまって助かりました。
たとえば
y'=xyとかだったらどうしよう？って話でした！


(PC版ではクリックすると原寸大が見れます、スマホ？ピンチアウト(スワイプみたいなあれ)でなんとかなるんじゃないすか(適当))
 

xの構成方法と、xに関するk1～k4の計算方法を図に記述しました。
y,zも同様です

この状態で、h=0.01、n=2000くらいまで計算させました。

それからですね、

こんな感じの3次行列を用意します。
「x軸回転」のすぐ下にdegreeの角度θ、そのさらに下にradの角度θを格納し、行列の意味は

こんな感じです。これをy軸回転、z軸回転で同様に繰り返し

このような感じに並べます。
一番下の3次行列はx,y,z、3軸回転行列の積で

このような意味合いです。
=MMULT(MMULT(x軸:x軸,y軸:y軸),z軸:z軸)
という演算をやっています。配列全体に染み渡らせたいときは、演算させたい範囲を囲ってｃｔｒｌ＋shift＋エンターです。
なお、sinの符号は、片方がもう片方の逆符号であればどっちがプラスでもいいです。どうせ逆回転になるだけなんで。


そうして出来上がった3D回転行列を、
元のx,y,zに作用させましょう。
=MMULT(x:z,$合成行列)

そうして出来上がった新しいx,y,z座標のxとyだけをグラフに反映させると

このようになります。

x→y→z軸の順番で、180度ずつ回転させてますので、全部の回転が済めば元の回転角に戻ります


この図が動く仕組みを説明しましょう

まず
now()-today()
を入力します。
関数now()は時刻込み今のの日付のシリアル値を出力します。
関数today()は、時刻抜きの今日の日付のシリアル値を出力します。

シリアル値というのは、1900年1月1日の0時0分0秒を「0」とし、24時間を「1」と定義した表現方法なので、
たとえば同日の朝6時というのは0.25という数値になります。

６：００と打って、書式をシリアル値から数値に戻すと、0.25になっているはずです。
また、1900/2/1と「イコールなしで」打てば、1月1日から数えて32日目なので、書式を数値に戻すと32になるはずです。
これらを踏まえて、日付に半角スペースを入れてから時刻を入力すると、たとえば
1900/2/1 6:00:00
でしたら、32.25という数値になります。


now()はたとえば今が2019/6/4の10:57だったら
2019/6/4 10:57のシリアル値を算出するので

43620.45758

ぐらいの数値になります。

実はこの数値、一番小さい桁を見ると、ミリ秒単位ぐらいで動いています。
F4を連打するとうごめいているのがわかるかと思います。

これを動画の媒介変数に用いることができます。

たとえばnow()の値を700万倍に拡大してグラフに反映させることができれば、目に見えて動きます。
が、整数部が5桁もあるので、700万倍したらオーバーフローしてしまいます。

これを解決する手っ取り早い方法として、(色んな方法がありますが)「関数today()で関数now()を引き算する」というのがあげられると思います。

必ず1未満の数値になりますし、引き算という演算を行っているため、Excelはシリアル値ではなく数値と勝手に解釈して出力してくれます。


今回は、180度の回転を3回繰り返したいので、180×3=540°の周期が必要です。
ですので、now()-today()を700万倍した数値を540で割り算した剰余(あまり)が役に立ってくれます。
剰余の関数はモジュロ演算なので、modという関数です。
A=mod(700万*(now()-today()),540)
という演算を行います。


x,y,z軸の回転角θ1、θ2、θ3のdegreeにそれぞれ
θ1=IF(A<=180,A,180)
θ2=IF(A<=180,0,IF(A<=180*2,A-180,180))
θ3=IF(A<=180*2,0,A-180*2)

と入力すると、上のように動くグラフが作れます




バタフライ効果という名前の由来そのものが
・「蝶が舞えば気象兵器が儲かる」
と
・ローレンツアトラクタのグラフが蝶のように見える
の2種類あるというのがなんとも、
2つの過去が1つの未来に収束したみたいで興味深いです

ｼｭﾀｯL・ψ・卍=3ｹﾞｰ