Метод Дорманда-Принса¶
Определение¶
Одноэтапный расчет в методе Дорманда-Принса выполняется следующим образом.
Тогда значение следующего шага \(y_{k+1}\) вычисляется как
Это вычисление методом Рунге-Кутты порядка 4. Мы должны знать, что мы не используем \(k_2\), хотя оно используется для вычисления \(k_3\) и так далее.
Далее мы вычислим значение следующего шага \(z_{k+1}\) методом Рунге-Кутты порядка 5 как
Рассчитываем разницу двух следующих значений \(\left|z_{k+1}-y_{k+1}\right|\).
Это считается ошибкой в \(y_{k+1}\). Мы вычисляем оптимальный интервал времени \(h_{opt}\) как,
где \(h\) в правой части - старый интервал времени. В практическом программировании этот новый \(h_{opt}\) будет использоваться на следующем этапе расчета, хотя автор считает, что его также следует использовать в текущих расчетах, когда он очень мал, например, наполовину или меньше.
Использование¶
Представьте, что мы хотим минимизировать следующее дифференциальное уравнение:
Тогда код будет выглядеть так:
// 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;
}