UnityでVisual C#用の常微分方程式ソルバーOpen Solving Library for ODEs(OSLO)を使う(2)ブルースカイ・カタストロフィを生じるGavrilov Shilnikov modelを計算してDormand＆PrinceのRK547Mで計算して玉を動かして軌跡を付ける。ぐるぐる回っていたと思ったら突然広がって戻る。

今回は千葉逸人さんが中二病用語と言っていたこれをやってみる。

 Python＋Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う（その11）　中二病用語、ブルースカイ・カタストロフィを生じるGavrilov Shilnikov modelを計算してGIFアニメに。

詳しくはスカラーペディアに。

Blue-sky catastrophe

この式を計算してます。

早速動画から。

ぐるぐる回っていたと思ったら突然広がって戻るのが面白い。スカラーペディアに出ていた図ともよく合ってる。

ソースコードはこちら。

using UnityEngine;

using System;

using System.Collections.Generic;

using Microsoft.Research.Oslo;

using System.Linq;

public class MathNET05 : MonoBehaviour

{

    List<float> x0 = new List<float>();

    List<float> x1 = new List<float>();

    List<float> x2 = new List<float>();

    int count;

    int n;

    void Start()

    {

        double t0 = 0.0;

        double tmax = 1000.0;

        n = 20000;

        double dt = tmax / n;

        var sol = Ode.RK547M(

            0,

            new Vector(0.012277918, -2.356078578, 0.018241293),

            (t, x) => Bluesky(t, x),

            new Options

            {

                AbsoluteTolerance = 1e-12,

                RelativeTolerance = 1e-12

            });

        var points = sol.SolveFromToStep(t0, tmax, dt).ToArray();

        n = points.Length;

        count = 0;

        foreach (var s in points)

        {

            x0.Add(Convert.ToSingle(s.X[0]));

            x1.Add(Convert.ToSingle(s.X[1]));

            x2.Add(Convert.ToSingle(s.X[2]));

        }

    }

    void Update()

    {

        if (count < n)

        {

            transform.position = new Vector3(x2[count], x1[count], x0[count]);

        }

        count += 5;

    }

    Vector Bluesky(double t, Vector x)

    {

        double[] xdot = new double[3];

        double myu = 0.456;

        double eps = 0.0357;

        xdot[0] = x[0] * (2.0 + myu - 10.0 * (x[0] * x[0] + x[1] * x[1]))

            + x[2] * x[2] + x[1] * x[1] + 2.0 * x[1];

        xdot[1] = -x[2] * x[2] * x[2] - (1 + x[1]) * (x[2] * x[2] + x[1] * x[1] + 2.0 * x[1])

            - 4.0 * x[0] + myu * x[1];

        xdot[2] = (1 + x[1]) * (x[2] * x[2]) + x[0] * x[0] - eps;

        return new Vector(xdot);

    }

}