使用 Math.Net Numerics 求解常微分方程组

Question

我正在尝试使用 Math.Net Numerics 求解常微分方程组。这不是我遇到问题的数学。我要解决的系统是：

x'(t)=x(t)+2y(t)+2t

y'(t)=3x(t)+2y(t)-4t

x(0)=-7/4, y(0)=55/8

我想使用四阶 Runge Kutta ODE 求解器来求解这个系统： https://numerics.mathdotnet.com/api/MathNet.Numerics.OdeSolvers/RungeKutta.htm

它传入的第一个参数是

Vector<T>y0

我遇到的问题是当我点击vector<T>的链接时，我没有找到基本信息。所以我很难理解这门课。你了解它以及如何使用它吗？

接下来的参数没问题，直到我来到最后一个：Func<double, Vector<double>, Vector<double>>，因为我不懂向量类我不明白如何制作这个。你知道如何创建这个 Func 吗？

我遇到的最后一件事是返回类型。它说它返回一个Vector'1[]。但这是什么？我找不到它的任何文档。

你能帮帮我吗？也许有人可以展示如何解方程？

Answer 1

MathNet.Numerics 类型文档似乎缺少有关如何创建向量的信息，但是文档kindly linked by Frank Modica: 对此进行了解释。

要创建Vector<int>，请编写如下内容：

var vector = Vector<int>.Build.Dense(new[] { 1, 2, 3 });

这使用了其中一种初始化函数，看起来有许多不同的类型。

Vector<T> 是一个generic type，这意味着它的成员被声明为具有特定类型，但这可以是任何满足其类型约束的类型。您可以使用尖括号显式声明您希望泛型类型具有哪种类型，或者让编译器推断它们，这就是上面的示例所做的。此示例将泛型类型实例化显式声明为 Vector<int>，但在某些情况下，编译器也可以从参数的类型中推断出类型。

您需要知道的另一件事是如何声明Func<T1, T2, TResult>。库文档也不会解释这一点，因为它是另一个核心 C# 概念。这称为anonymous function。

要将匿名委托传递给求解器，请编写如下内容：

var result = SecondOrder(y0, start, end, N, (someValue, someVector) => { /* some operation on someValue and someVector which returns a new Vector<T>*/ });

我不熟悉 Runge-Kutta，所以我不知道您的匿名方法的实现应该是什么样子，但我已经用注释对其进行了注释。

表达式：

(someValue, someVector) => { return answer; } 是一个 Lambda 表达式，编译器将其理解为与预期的 Func<double, Vector<double>, Vector<double>> 匹配 - 它声明了一个函数，该函数接受一个 double 和一个向量，并返回一个向量。

最后，Vector'1[] 语法令人困惑是可以理解的——这似乎是文档格式的问题。这就是编译器在内部表示匿名类型的方式，但不应该真正暴露给用户。它与Vector<T>[] 的含义相同，或者泛型向量数组。

Answer 2

我想我解决了。感谢 Tom W 的帮助。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.OdeSolvers;

namespace Sys
{
    class Program
    {
        static void Main(string[] args)
        {
            int N = 1000000;
            Vector<double> y0 = Vector<double>.Build.Dense(new[] { -7.0 / 4.0, 55.0 / 8.0 });
            Func<double, Vector<double>, Vector<double>> der = DerivativeMaker();

            Vector<double>[] res = RungeKutta.FourthOrder(y0, 0, 10, N, der);

            double[] x = new double[N];
            double[] y = new double[N];
            for (int i=0; i <N; i++)
            {
                double[] temp = res[i].ToArray();
                x[i] = temp[0];
                y[i] = temp[1];
            }

            //Test
            Console.WriteLine(y[N / 10]); // gives 164,537981852489
            Console.WriteLine(Math.Exp(-1) + 3 * Math.Exp(4) - 5.0 / 2 + 23.0 / 8); //gives 164,537329540604, which is y(1)


            Console.ReadKey();

        }

        static Func<double, Vector<double>, Vector<double>> DerivativeMaker()
        {
            return (t, Z) =>
            {
                double[] A = Z.ToArray();
                double x = A[0];
                double y = A[1];

                return Vector<double>.Build.Dense(new[] { x + 2 * y + 2 * t, 3 * x + 2 * y - 4 * t });

            };
        }
    }
}