Метод Дорманда-Принса

Определение

Одноэтапный расчет в методе Дорманда-Принса выполняется следующим образом.

\begin{equation} \begin{aligned} k_{1} &=h f\left(t_{k}, y_{k}\right), \\ k_{2} &=h f\left(t_{k}+\frac{1}{5} h, y_{k}+\frac{1}{5} k_{1}\right), \\ k_{3} &=h f\left(t_{k}+\frac{3}{10} h, y_{k}+\frac{3}{40} k_{1}+\frac{9}{40} k_{2}\right), \\ k_{4} &=h f\left(t_{k}+\frac{4}{5} h, y_{k}+\frac{44}{45} k_{1}-\frac{56}{15} k_{2}+\frac{32}{9} k_{3}\right), \\ k_{5} &=h f\left(t_{k}+\frac{8}{9} h, y_{k}+\frac{19372}{6561} k_{1}-\frac{25360}{2187} k_{2}+\frac{64448}{6561} k_{3}-\frac{212}{729} k_{4}\right), \\ k_{6} &=h f\left(t_{k}+h, y_{k}+\frac{9017}{3168} k_{1}-\frac{355}{33} k_{2}-\frac{46732}{5247} k_{3}+\frac{49}{176} k_{4}-\frac{5103}{18656} k_{5}\right), \\ k_{7} &=h f\left(t_{k}+h, y_{k}+\frac{35}{384} k_{1}+\frac{500}{1113} k_{3}+\frac{125}{192} k_{4}-\frac{2187}{6784} k_{5}+\frac{11}{84} k_{6}\right). \end{aligned} \end{equation}

Тогда значение следующего шага \(y_{k+1}\) вычисляется как

\begin{equation} y_{k+1}=y_{k}+\frac{35}{384} k_{1}+\frac{500}{1113} k_{3}+\frac{125}{192} k_{4}-\frac{2187}{6784} k_{5}+\frac{11}{84} k_{6}. \end{equation}

Это вычисление методом Рунге-Кутты порядка 4. Мы должны знать, что мы не используем \(k_2\), хотя оно используется для вычисления \(k_3\) и так далее.

Далее мы вычислим значение следующего шага \(z_{k+1}\) методом Рунге-Кутты порядка 5 как

\begin{equation} z_{k+1}=y_{k}+\frac{5179}{57600} k_{1}+\frac{7571}{16695} k_{3}+\frac{393}{640} k_{4}-\frac{92097}{339200} k_{5}+\frac{187}{2100} k_{6}+\frac{1}{40} k_{7} \end{equation}

Рассчитываем разницу двух следующих значений \(\left|z_{k+1}-y_{k+1}\right|\).

\begin{equation} \left|z_{k+1}-y_{k+1}\right|=\left|\frac{71}{57600} k_{1}-\frac{71}{16695} k_{3}+\frac{71}{1920} k_{4}-\frac{17253}{339200} k_{5}+\frac{22}{525} k_{6}-\frac{1}{40} k_{7}\right| \end{equation}

Это считается ошибкой в \(y_{k+1}\). Мы вычисляем оптимальный интервал времени \(h_{opt}\) как,

\begin{equation} s =\left(\frac{\varepsilon h}{2\left|z_{k+1}-y_{k+1}\right|}\right)^{\frac{1}{5}}, \\ h_{o p t} =s h, \end{equation}

где \(h\) в правой части - старый интервал времени. В практическом программировании этот новый \(h_{opt}\) будет использоваться на следующем этапе расчета, хотя автор считает, что его также следует использовать в текущих расчетах, когда он очень мал, например, наполовину или меньше.

Использование

Представьте, что мы хотим минимизировать следующее дифференциальное уравнение:

\begin{equation} y' = 3 \frac{y}{x} + x^3 + x, y(1) = 3 \end{equation}

Тогда код будет выглядеть так:

// example_dorpi.cpp

#include <iostream>
#include "../src/numerary.hpp" // Numerary library

using namespace std;
using namespace numerary;

/* Equation to solve */
double equation(double x, double y) {
    return 3.0*y/x + x*x*x + x;
}

/* The main function */
int main() {

    double *y = new double[2];
    double x0, x, h;
    short int is_solved;

    // Initial point
    x0 = 1; y[0] = 3;

    // Point where we want calculate y(x)
    x = 2.0;

    // Step size
    h = 0.01;

    is_solved = Numerary::ordinary_differential_equations(equation, y, x0, h, x, "dorpi_4_5");

    if (is_solved == 0) {
        cout << "Solved!" << endl;
        cout << "y(" << x << ") = " << y[1] << endl;
    } else {
        cout << "Method not allowed!" << endl;
    }

    delete[] y;

    return 0;
}