Язык программирования C#

Решение задачи Коши при помощи методов Рунге-Кутты

Рассматривается задача Коши вида

\[ \left\{ \begin{array} \mathbf{u}'= \mathbf{f}(x, \mathbf{u}),\\ \mathbf{u}(x_0)=\mathbf{u}_0, \end{array} \right. \]

где \( \mathbf{u} \) - вектор неизвестных функций, \( x \) - переменная, \( \mathbf{f} \) - набор вектор-функций нескольких переменных.

где введено обозначение

\[ \begin{array}{l} \mathbf{k}_1 = h\mathbf{f}(x,u(x))\\ \mathbf{k}_2 = h\mathbf{f}\left(x+\alpha_2h,u(x)+\beta_{21}\mathbf{k}_1\right)\\ .....................................\\ \mathbf{k}_q = h\mathbf{f}\left(x+\alpha_qh,u(x)+\beta_{q1}\mathbf{k}_1+...+\beta_{q,q-1}\mathbf{k}_{q-1}\right) = h\mathbf{f}\left(x+\alpha_qh,u(x)+\sum\limits_{i=1}^{q-1}\beta_{qi}\mathbf{k}_i\right)\\ \end{array} \]

Конкретный метод определяется числом \( q \) и коэффициентами \( \alpha_i \) и \( \beta_{ij} \). Эти коэффициенты часто упорядочивают в таблицу (называемую таблицей Батчера):

\[ \begin{array}{c|ccccc}  0 &&&&&\\ \alpha_2 & \beta_{21} &&&&\\ \alpha_3 & \beta_{31} & \beta_{32} &&&\\ \vdots & \vdots & \vdots& \ddots&&\\ \alpha_q & \beta_{q1} & \beta_{q2}& \dots & \beta_{q,q-1}&\\ \hline & p_1 & p_2 & \dots & p_{q-1} & p_q \end{array} \]

Вложенные формулы Рунге-Кутты

Используются следующие формулы

\[ u(x+h)=u(x)+\sum\limits_{i=1}^q p_ik_i \]

\[ \tilde u(x+h)=u(x)+\sum\limits_{i=1}^q \tilde p_ik_i \]

Формула для \( u(x+h) \) имеет погрешность порядка \( p \), Формула для \( \tilde u(x+h) \) имеет погрешность \( s = p\pm 1 \). Для практической оценки погрешности получается выражение

\[ r(x+h) = u(x+h) - \tilde u(x+h)\approx\sum\limits_{i=1}^q \left(p_i-\tilde p_i\right)k_i \]

Таблица Батчера приобретает вид:

\[ \begin{array}{c|ccccc}  0 &&&&&\\ \alpha_2 & \beta_{21} &&&&\\ \alpha_3 & \beta_{31} & \beta_{32} &&&\\ \vdots & \vdots & \vdots& \ddots&&\\ \alpha_q & \beta_{q1} & \beta_{q2}& \dots & \beta_{q,q-1}&\\ \hline & p_1 & p_2 & \dots & p_{q-1} & p_q \\ & \tilde p_1 & \tilde p_2 & \dots & \tilde p_{q-1} & \tilde p_q \end{array} \]

Решение задачи Коши средствами языка C#

1. Для решения задачи создаём новый проект консольного приложения под названием RungeKutta и добавляем в проект новый обобщенный интерфейс (в отдельном файле).

       /// <summary>
       /// Интерфейс для решения дифференциальных уравнений
       /// </summary>
       /// <typeparam name="T">тип переменной</typeparam>
       /// <typeparam name="U">тип значения функции</typeparam>
       interface IOdEsolver<T, U>
       {
           /// <summary>
           /// правые части дифференциальных уравнений
           /// </summary>
           List<Func<T, List<U>, U>> Equations { get; } // List<T>, Func<in T1, in T2, out TResult> 
   
           /// <summary>
           /// погрешность вычислений
           /// </summary>
           double Epsilon { get; }
   
           /// <summary>
           /// метод для решения задачи Коши
           /// </summary>
           /// <param name="a">точка, в которой заданы начальные условия</param>
           /// <param name="b">точка, в которой ищется решение</param>
           /// <param name="initialConditions">вектор начальных условий</param>
           /// <returns>решение задачи Коши в точке x=b</returns>
           List<T> Solve(T a, T b, List<T> initialConditions);
       }
   

2. Добавляем в проект статический класс ButcherTableau, содержащий таблицы Батчера для нескольких методов из семейства методов Рунге-Кутты: метода Хойна, метода Богацки-Шампайна, метода Кеша-Карпа, метода Рунге-Кутты-Фельдберга и двух методов Дорманда-Принса разных порядков точности.
3. Добавляем в проект класс, реализующий интерфейс IOdEsolver

   public class OdeSolver : IOdEsolver<double, double>
   

4. Добавляем в новый класс конструктор

   public OdeSolver(List<Func<double, List<double>, double>> funcs, double[][] butcherTableau, double epsilon = 0.1e-6)
   

5. Реализуем метод Solve. Алгоритм его работы следующий
   1. Вычисляем новую переменную h = b - a
   2. В цикле, пока a<b (используем условие Math.Abs(b - a) > double.Epsilon)
   1. Вычисляем \( u(a+h) \) и \( \tilde u(a+h) \).
   2. Вычисляем погрешность \( r(a+h) \)
   3. Если модуль погрешность превышает Epsilon, делим h пополам и снова вычисляем \( u(a+h) \) и \( \tilde u(a+h) \).
   4. a += h;
   5. Вектор начальных условий приравниваем \( \tilde u(a+h) \)
   3. Возвращаем \( \tilde u(a+h) \).
6. Создаём модульные тесты. Код тестирующего метода должен быть примерно таким

           [TestMethod()]
           public void SolveTest()
           {
               var equation = new List<Func<double, List<double>, double>> { (x, y) => y[0] };
               OdeSolver odeSolver = new OdeSolver(equation, ButcherTableau.RungeKuttaFeldberg, 0.1e-6);
               var solution = odeSolver.Solve(0, 1, new List<double> { 1 });
               Assert.IsTrue(Math.Abs(solution[0] - Math.E) < 0.1e-5);
   
               var equations = new List<Func<double, List<double>, double>> {(x, y) => y[1], (x, y) => -y[0]};
               OdeSolver solver = new OdeSolver(equations, ButcherTableau.RungeKuttaFeldberg, 0.1e-6);
               solution = solver.Solve(0, 1, new List<double> {0, 1});
               Assert.IsTrue(Math.Abs(solution[0] - Math.Sin(1)) < 0.1e-5);
               Assert.IsTrue(Math.Abs(solution[1] - Math.Cos(1)) < 0.1e-5);
               
               // добавьте свои тесты для каждого из методов 
           }
   

7. Продемонстрируйте использование класса OdeSolver в основной программе