Постановка задачи Коши

        Метод Рунге-Кутта относится к методам численного решения обычных дифференциальный уравнений первого порядка. Этот численный метод является одним из точных методов численного решения ОДУ и систем ОДУ, и не очень сложных в реализации, поэтому этот медот получил широкое распространение. Задача Коши для ОДУ первого порядка ставится следующим образом:

dy/dx=F(x,y)

        Как и для любых других чесленных методов решения дифференциальных уравнений, для решения этого уравнения требуется задать начальные условия: x0, y0.

Математическое решение задачи

        Решением поставленной задачи будет ряд точек на плоскости (x,y). Обозначим шаг вычислений как h. Вектор точек по оси x обозначим x[i], i=0...N, при этом значения этого вектора будут определяться следующим образом: x[i]=x[0]+h*i. Вектор точек по оси y обозначим как y[i]. Тогда значения ветора y будут определяться следующим образом:

y[i+1]=y[i]+delta;
delta=(K1+2*K2+2*K3+K4)/6;
K1=h*F(x[i],y[i]);
K2=h*F(x[i]+h/2,y[i]+K1/2);
K3=h*F(x[i]+h/2,y[i]+K2/2);
K4=h*F(x[i]+h,y[i]+K3);

        Таким образом будет получено численное решение исходного дифференциального уравнения на интервале [a,b] с заданными начальными условиями и шагом.

Реализация метода на C++

        Оформим метод в виде функции на C++. Фактически вся "начинка" уже написана выше. Однако не будем искать "лёгких" путей, и напишем так, чтобы получившийся код был удобен в использовании и работал быстро. Для начала определим новый тип данных:

typedef std::pair<double, double> pair_of_doubles;

        Этот тип данных представляет собой пару значений типа double, которая будет хранить значение по x в поле first и значение типа y в поле second. Далее определим вспомогательную функцию RungeKuttStep, которая будет делать шаг по методу Рунге-Кутта (по описанным выше формулам). Эта функция, с точки зрения пользоватля, может пригодится для нахождения оптимального шага и оценки погрешности по правилу Рунге.

double RungeKuttStep(double x, double y, double h, double (*F)(double x, double y))
{
double K1 = h*F(x ,y );
double K2 = h*F(x+h/2.0,y+K1/2.0);
double K3 = h*F(x+h/2.0,y+K2/2.0);
double K4 = h*F(x+h ,y+K3 );
return y+(K1+2.0*K2+2.0*K3+K4)/6.0;
}

        Давайте разберёмся, что функция принимает и что возвращает (надеюсь то, что она делает Вам уже понятно - см. выше). Функция принимает текущее значение x, текущее значение y, значение шага h, и, самое главное, указатель на функцию вида double F(double x, double y). Функция F будет вычислять значение функции при соответствующих x и y. Функция RungeKuttStep вернёт полученное значение y при x+h.
        Далее определим основную функцию RungeKutt, которая будет вычислять набор значений (x,y) в заданном интервале. Для этого определим её следующим образом:

bool RungeKutt(double x0, double y0, double xn, double h, double (*F)(double x, double y), std::vector<pair_of_doubles> &result)
{
result.clear();          // Очищаем вектор результатов, вдруг пользователь не позаботился об этом
int n = ceil((xn-x0)/h); // Считаем необходимое количество шагов
result.reserve(n+1);     // А резервиреум на одно больше, т.к. понадобится место для (x0,y0)
double x = x0;           // Текущие значение x и y равны
double y = y0;           // соответственно x0 и y0
result.push_back(std::pair<double,double>(x,y)); // Сохраняем первую точку
for(int i=0;i<n;++i)
 { // Последовательно выполняем полученное количество шагов
 y = RungeKuttStep(x,y,h,F);
 x+=h; // Не забывая "продвигать" x и сохранять результаты
 result.push_back(std::pair<double,double>(x,y));
 }
return true;
}

        В этой функции входные параметры похожи на параметры функции RungeKuttStep: x0,y0 - начальные условия, h - шаг, F - указатель на функцию. Добавились лишь параметры xn - правая граница интервала (исследование проводится на интервале [x0,xn]) и result - вектор результатов. после вызова функции Вам остаётся лишь перенести полученные точки на график.

Материалы к статье