UnityでVisual C#用の常微分方程式ソルバーOpen Solving Library for ODEs(OSLO)を使ってピタゴラスの三体問題を適応型ルンゲクッタ法、Dormand＆PrinceのRK547Mで計算して玉を動かして軌跡を付ける。カメラも回転させる。ついでに背景もつけた。

UnityでMath.NET Numericsという数値計算ライブラリを使うシリーズをやっているが、ルンゲクッタ法が4段4次のものしかないのでピタゴラスの三体問題のような精度が問題になるようなものはできず、Pythonで計算したものをCSVファイルにして読み込むということしかできなかった。

しかし、MicrosoftがOpen Solving Library for ODEs(OSLO)を出していることを知ったのでこれならできるな、とやってみた。

AssetsフォルダにMicrosoft.Research.Oslo.dllを置くと自動的に参照してくれた。

ソースコードは以下の通り。

using UnityEngine;

using System;

using System.Collections.Generic;

using Microsoft.Research.Oslo;

using System.Linq;

public class CsvReader : MonoBehaviour

{

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

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

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

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

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

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

    int n;

    int count;

    GameObject Ball1;

    GameObject Ball2;

    GameObject Ball3;

    // Start is called once before the first execution of Update after the MonoBehaviour is created

    void Start()

    {

        Ball1 = GameObject.Find("Ball1");

        Ball2 = GameObject.Find("Ball2");

        Ball3 = GameObject.Find("Ball3");

        var sol = Ode.RK547M(

                                0,

                                new Vector(1, 0, 3, 0, -2, 0, -1, 0, 1, 0, -1, 0),

                                (t, x) => threebody(t, x),

                                new Options

                                {

                                    AbsoluteTolerance = 1e-12,

                                    RelativeTolerance = 1e-12

                                });

        var points = sol.SolveFromToStep(0, 100, 0.01).ToArray();

        n = points.Length;

        count = 0;

        foreach (var s in points)

        {

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

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

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

            y2.Add(Convert.ToSingle(s.X[6]));

            x3.Add(Convert.ToSingle(s.X[8]));

            y3.Add(Convert.ToSingle(s.X[10]));

        }

    }

    Vector threebody(double t, Vector x)

    {

        var f = new double[12];

        var m1 = 3.0;

        var m2 = 4.0;

        var m3 = 5.0;

        var qx1 = x[0];

        var vx1 = x[1];

        var qy1 = x[2];

        var vy1 = x[3];

        var qx2 = x[4];

        var vx2 = x[5];

        var qy2 = x[6];

        var vy2 = x[7];

        var qx3 = x[8];

        var vx3 = x[9];

        var qy3 = x[10];

        var vy3 = x[11];

        var Ra = Math.Pow(Math.Sqrt(Math.Pow((qx2 - qx1), 2.0) + Math.Pow((qy2 - qy1), 2)), 3);

        var Rb = Math.Pow(Math.Sqrt(Math.Pow((qx3 - qx1), 2.0) + Math.Pow((qy3 - qy1), 2)), 3);

        f[0] = vx1;

        f[1] = m2 * (qx2 - qx1) / Ra + m3 * (qx3 - qx1) / Rb;

        f[2] = vy1;

        f[3] = m2 * (qy2 - qy1) / Ra + m3 * (qy3 - qy1) / Rb;

        Ra = Math.Pow(Math.Sqrt(Math.Pow((qx1 - qx2), 2.0) + Math.Pow((qy1 - qy2), 2)), 3);

        Rb = Math.Pow(Math.Sqrt(Math.Pow((qx3 - qx2), 2.0) + Math.Pow((qy3 - qy2), 2)), 3);

        f[4] = vx2;

        f[5] = m1 * (qx1 - qx2) / Ra + m3 * (qx3 - qx2) / Rb;

        f[6] = vy2;

        f[7] = m1 * (qy1 - qy2) / Ra + m3 * (qy3 - qy2) / Rb;

        Ra = Math.Pow(Math.Sqrt(Math.Pow((qx1 - qx3), 2.0) + Math.Pow((qy1 - qy3), 2)), 3);

        Rb = Math.Pow(Math.Sqrt(Math.Pow((qx2 - qx3), 2.0) + Math.Pow((qy2 - qy3), 2)), 3); ;

        f[8] = vx3;

        f[9] = m1 * (qx1 - qx3) / Ra + m2 * (qx2 - qx3) / Rb;

        f[10] = vy3;

        f[11] = m1 * (qy1 - qy3) / Ra + m2 * (qy2 - qy3) / Rb;

        return f;

    }

    // Update is called once per frame

    void Update()

    {

        if (count > n)

        {

            count = 0;

        }

        Ball1.transform.position = new Vector3(x1[count], y1[count], 0f);

        Ball2.transform.position = new Vector3(x2[count], y2[count], 0f);

        Ball3.transform.position = new Vector3(x3[count], y3[count], 0f);

        count++;

    }

}

結果をカメラを回しながら動画にしたもの。せっかくなんで背景に空もつけた。

ちゃんと計算してそうで一安心。