Visual C# (C_sharp)の数学ライブラリ Math.NET Numericsを使う(5) 常微分方程式の数値解法、4段4次のルンゲクッタ法がRungeKutta.FourthOrderの一文でできる。ローレンツ方程式を例としてやってみる。

C#のMath.NET Numericsを使った数値解析シリーズも5回目。さて今回は常微分方程式の数値解法。2次と4次のルンゲクッタ法とアダムスバッシュフォース法が使える。 まずは最初にusing MathNet.Numerics.OdeSolvers;はつけておく。で基本、RungeKutta.FourthOrderの1文で計算できる。Vectorの配列に結果は入るのと、関数もFunc(double, Vector,Vector)で与えることだけ注意。

例題はまあLorenz方程式がいいでしょう。

\begin{align} &\frac{dx}{dt} = \sigma (y - x) \\\ &\frac{dy}{dt} = -xz + rx -y \\\ &\frac{dz}{dt} = xy - bz \\\ \end{align}

プログラムはこんな感じで、

結果はこちら：

ソースコードの詳細はこちら：

using System;
using System.Windows.Forms;
using System.Windows.Forms.DataVisualization.Charting;
using MathNet.Numerics;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.OdeSolvers;
using Series = System.Windows.Forms.DataVisualization.Charting.Series;

namespace RungeKuttaTest01
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
            chart1.Series.Clear();
            chart1.Titles.Clear();
            chart1.Legends.Clear();
            chart1.ChartAreas.Clear();
        }

        private void button1_Click(object sender, EventArgs e)
        {
            Title title = new Title("Lorenz Equation", Docking.Top);
            chart1.Titles.Add(title);
            Legend legend = new Legend();
            chart1.Legends.Add(legend);

            Series series1 = new Series();
            series1.ChartType = SeriesChartType.Line;
            series1.BorderWidth = 1;
            series1.LegendText = "Y";
            chart1.Series.Add(series1);

            Series series2 = new Series();
            series2.ChartType = SeriesChartType.Line;
            series2.BorderWidth = 1;
            series2.LegendText = "Z";
            chart1.Series.Add(series2);
            chart1.ChartAreas.Add("");
            Axis axisX = new Axis();
            axisX.Title = "X";
            axisX.Minimum = -20;
            axisX.Maximum = 25;
            axisX.Interval = 5;
            chart1.ChartAreas[0].AxisX = axisX;

            Axis axisY = new Axis();
            axisY.Title = "Y,Z";
            axisY.Minimum = -30;
            axisY.Maximum = 60;
            axisY.Interval = 10;
            chart1.ChartAreas[0].AxisY = axisY;

            double dt = 0.01;
            double tmax = 200.0;
            int n = Convert.ToInt32(tmax / dt);

            var t = Generate.LinearSpaced(n, 0.0, tmax);
            var x0 = Vector.Build.DenseOfArray(new double[] { 1.0, 1.0, 1.0 });
            var x = RungeKutta.FourthOrder(x0, 0, tmax, n, Func);


            for (int i = 0; i < n; i++)
            {
                series1.Points.AddXY(x[i][0], x[i][1]);
                series2.Points.AddXY(x[i][0], x[i][2]);
            }
        }

        public static Vector Func(double t, Vector x)
        {
            double s = 10.0, r = 28.0, b = 8.0/3.0;

            double x_dot = s * (x[1] - x[0]);
            double y_dot = r * x[0] - x[1] - x[0] * x[2];
            double z_dot = x[0] * x[1] - b * x[2];

            return Vector.Build.DenseOfArray(new double[] { x_dot, y_dot, z_dot });
        }
    }
}

過去のもの：

Visual C# (C_sharp)の数学ライブラリ Math.NET Numericsを使う(1) 複素行列を定義して一次方程式や逆行列、行列式などを計算する。

Visual C# (C_sharp)の数学ライブラリ Math.NET Numericsを使う(2) 補間を行う（Interpolate) リニア、３次スプライン、有理関数などいろいろ使える。

Visual C# (C_sharp)の数学ライブラリ Math.NET Numericsを使う(3) 高速フーリエ変換（FFT）を実行する。FourierOptionsにMatlabとNumerical Recipesがあるのが意外。

Visual C# (C_sharp)の数学ライブラリ Math.NET Numericsを使う(4) 多項式フィッティングをして、Array.ConvertAllで一括でフィッティングデータを得る。

Visual C#でMath.NET Numericsを使うシリーズ全編はこちらから

高周波＆数値計算関係記事リンク集