Visual C#用のMicrosoftの常微分方程式ソルバー、Open Solving Library for ODEs（OSLO）を使ってみる(1) 適応型のルンゲクッタ、Dormand＆PrinceのRK547MとGear の後退差分式(BDF)が使える。ピタゴラスの三体問題を計算してみた。

Visual C#で数値計算する時はMath.NET Numericsを使うのが便利だが、常微分方程式が4段4次のルンゲクッタ法しかないので硬い問題などがほぼ無理。そこでもうちょっといろいろあるのはないかなと思って探すとこれがあった。

Open Solving Library for ODEs

適応型のルンゲクッタ、Dormand＆PrinceのRK547MとGear の後退差分式(BDF)が使える。まずはダウンロードして、、、

あれ？dllがない？自分でビルドしろということか。早速やってみると（古いバージョンなので）いろいろ怒られながらとりあえずできた。

使い方はとても簡単なので、早速やってみる。4段4次のルンゲクッタで計算するのがとても難しいピタゴラスの三体問題にしよう。

 Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う（その5）　ピタゴラスの三体問題を計算する。rtolとatolを設定しないと無茶苦茶になる。

Visual C#のコードはこんな感じで。RK547Mを使った。

using Microsoft.Research.Oslo;

namespace PlotRungeKutta02

{

    public partial class Form1 : Form

    {

        public Form1()

        {

            InitializeComponent();

            var axis = formsPlot1.Plot.Axes;

            axis.Bottom.Label.Text = "x";

            axis.Bottom.TickLabelStyle.FontSize = 21;

            axis.Left.Label.Text = "y";

            axis.Left.TickLabelStyle.FontSize = 21;

            axis.Title.Label.Text = "3 body probrem";

            axis.Title.Label.FontSize = 32;

            axis.SetLimits(-6, 6, -6, 6);

            formsPlot1.Refresh();

        }

        private void button1_Click(object sender, EventArgs e)

        {

            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();

            var x1 = new double[points.Length];

            var y1 = new double[points.Length];

            var x2 = new double[points.Length];

            var y2 = new double[points.Length];

            var x3 = new double[points.Length];

            var y3 = new double[points.Length];

            int i = 0;

            foreach (var s in points)

            {

                x1[i] = s.X[0];

                y1[i] = s.X[2];

                x2[i] = s.X[4];

                y2[i] = s.X[6];

                x3[i] = s.X[8];

                y3[i] = s.X[10];

                i++;

            }

            var p1 = formsPlot1.Plot.Add.ScatterLine(x1, y1);

            p1.LineWidth = 5;

            var p2 = formsPlot1.Plot.Add.ScatterLine(x2, y2);

            p2.LineWidth = 5;

            var p3 = formsPlot1.Plot.Add.ScatterLine(x3, y3);

            p3.LineWidth = 5;

            formsPlot1.Refresh();

        }

        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;

        }

    }

}

結果はこちら。

ちゃんと計算できてるっぽい。Ode.RK547Mの代わりにOde.GearBDFにしても同様に計算できた。これは使えそう。