Алгоритм реализован в виде нескольких программ на языке С, которые вы можете скачать в разделе программы.

Опишем сначала, что такое дифференциальные уравнения на примере из механики и введем несколько обозначений.

Представьте себе, что вы бросили камень (можно бросать и гранату, но это негуманно) по направлению вверх и вперед. Из житейского опыта известно, что он полетит вперед по гладкой дуге, сначала вверх, потом вниз. На уроках физики вы может быть решали эту задачу аналитически и узнали, что камень летит по параболе. Естественно, для простоты будем считать, что воздух не оказывает никакого сопротивления камню.

На камень в полете действует только одна сила - сила тяжести. Если масса камня m, то сила тяжести, равна F=mg, где g-ускорение свободного падения, вблизи Земли g=9.81 (м с²)^-1

Введем координаты камня. Пусть x - расстояние по горизонтали от места с которого кидали камень, пусть y - расстояние от поверхности Земли до камня. Условно примем, что в начальный момент времени x=0, y=0. Если мы хотим узнать поведение камня в любой момент времени, естественно, что эти координаты зависят от времени: x(t), y(t). Пусть v_x(t), v_y(t) - скорость камня в горизонтальном и вертикальном направлении, соответственно. Пусть a_x(t), a_y(t) - ускорение камня в горизонтальном и вертикальном направлении соответственно.

Воспользуемся вторым законом Ньютона, который связывает ускорение движущегося тела и силы, действующие на него:

F = ma

Здесь F - сумма всех сил, действующих на тело, a - полное ускорение. Обе эти величины векторные. т.е. имеют направление.

Это же уравнение можно записать в проекциях на оси x, y.

F_x = m a_x - сила действующая в горизонтальном направлении.
F_y = m a_y - сила действующая в вертикальном направлении.

Сила тяжести действует на тело в вертикальном направлении и направлена вниз. То есть

F_x = m a_x = 0
F_y = m a_y = - m g

Как известно из курса кинематики скорость, ускорение и координата тела связаны друг с другом. Скорость равна производной координаты по времени. Ускорение равно производной скорости по времени, т.е. ускорение равно второй производной координаты по времени. Далее штрих (') означает производную по времени.

v_x = x'
v_y = y'
a_x = x''
a_y = y''

Заменим в уравнениях Ньютона ускорения на скорости и запишем все уравнения в виде системы, где в левой части стоят производные, а в правой части сами переменные.

m v'_x = 0
m v'_y = - m g
x' = v_x
y' = v_y

Поделим первое и второе уравнение на m:

x' = v_x
y' = v_y
v'_x = 0
v'_y = - g

Полученная система является системой линейных уравнений первой степени. Неизвестными в ней являются x(t), y(t), v_x(t), v_y(t). Для того, чтобы решение этой системы существовало и было единственным необходимо и достаточно, чтобы были заданы значения неизвестных в начальный момент времени t=0. Действительно, мы их чаще всего знаем. В нашей задаче мы так ввели оси координат, что x(0)=0, y(0)=0. Пусть мы знаем и проекции начальной скорости тела v_x(0), v_y(0). Чаще всего задают не проекции, а полную скорость тела v и угол к горизонту, под которым тело бросали.

Введем четырехмерный вектор U(x, y, v_x, v_y). Это может показаться странным, что для такой простой задачи нам понадобились какие-то четырехмерные вектора. Не пугайтесь. Четырехмерных векторов в природе не бывает. Это лишь абстракция, которая поможет нам проще записать решение этой задачи и программу, которая этот метод реализует. Итак, координатами этого вектора являются координаты камня и проекции его скорости. Нашу систему уравнений можно записать в виде:

U' = F(U,t)

где F(U,t) = (v_x, v_y, 0, -g)

Перейдем к описанию метода решения таких задач. Метода Рунге-Кутта.

Это численный метод. Суть метода состоит в том, чтобы разбить интервал времени, в течение которого нам надо узнать координаты и скорости на некоторое большое количество отрезков времени dt.

Пусть мы знаем значение вектора U в некоторый момент времени t и хотим узнать его в момент времени t+dt.

Методом Эйлера называется метод в котором
U(t+dt) = U(t) + F(U,t)dt.

Этот метод действительно работает. Посмотрите покомпонентно расписанную эту же запись:

x(t+dt) = x(t) + v_xdt
y(t+dt) = y(t) + v_ydt
v_x(t+dt) = v_x(t) + a_xdt
v_y(t+dt) = v_y(t) + a_ydt

Эти равенства являются приближенными и следуют из определения производной, Чем меньше dt, тем точнее эти равенства.

Методом Рунге-Кутта 4-го порядка называется метод в котором
U(t+dt) = U(t) + (k₁ + 2k₂ + 2k₃ + k₄)/6.
где

k₁ = F(U,t)dt
k₂ = F(U+k₁/2, t+dt/2)dt
k₃ = F(U+k₂/2, t+dt/2)dt
k₄ = F(U+k₃, t+dt)dt

Понимаю, что это довольно абстрактное выражение, непонятное интуитивно, но можете считать, что это тот же метод Эйлера, только уточненный. Более подробно о методе Рунге-Кутта (с доказательствами) можно почитать в книгах по численным методам (например Калиткин Н.Н., Численные методы).

Напишем часть кода на С++, который позволит нам реализовать этот метод:

class vector{
   public:
   float x, y, vx, vy;
   vector(){
       x = 0;
       y = 0;
       vx = 0;
       vy = 0;
   }
   vector(float ax, float ay, float avx, float avy){
       x = ax;
       y = ay;
       vx = avx;
       vy = avy;
   }
   void out(){
       printf("x=%f\t", x);
       printf("y=%f\t", y);
       printf("vx=%f\t", vx);
       printf("vy=%f\t", vy);
   }
};
vector operator +(vector a, vector b){
   vector c;
   c.x = a.x + b.x;
   c.y = a.y + b.y;
   c.vx = a.vx + b.vx;
   c.vy = a.vy + b.vy;
   return(c);
}
vector operator -(vector a, vector b){
   vector c;
   c.x = a.x - b.x;
   c.y = a.y - b.y;
   c.vx = a.vx - b.vx;
   c.vy = a.vy - b.vy;
   return(c);
}
vector operator *(float a, vector b){
   vector c;
   c.x = a * b.x;
   c.y = a * b.y;
   c.vx = a * b.vx;
   c.vy = a * b.vy;
   return(c);
}
vector operator *(vector b, float a){
   vector c;
   c.x = a * b.x;
   c.y = a * b.y;
   c.vx = a * b.vx;
   c.vy = a * b.vy;
   return(c);
}

В данном участке кода мы ввели класс под названием vector. У него есть 4 компонента: x, y, vx, vy. Далее мы ввели операции сложения и вычитания таких векторов (покомпонентное сложение или вычитание) а также операции умножения вектора на число и числа на вектор (которые совпадают с точностью до перестановки множителей местами).

Напишем функцию F(U,t) для нашей задачи с камнем:

vector F(vector U, float t){
   vector res;
   float g = 9.81;
   res.x = U.vx;
   res.y = U.vy;
   res.vx = 0;
   res.vy = - g;
   return(res);
}

Дальше нам осталось только написать главную часть программы, которая будет использовать этот код:

#include <stdio.h>
void main(){
   float t=0, dt = 0.05;
   float vx = 1, vy = 2;
   float x=0, y=0;
   vector U(x, y, vx, vy);
   vector k1, k2, k3, k4;
   while(U.y >= 0){
       k1 = F(U, t)*dt;
       k2 = F(U + 0.5*k1, t+0.5*dt)*dt;
       k3 = F(U + 0.5*k2, t+0.5*dt)*dt;
       k4 = F(U + k3, t+dt)*dt;
       U = U + 1.0/6.0 * (k1 + 2*k2 + 2*k3 + k4);
       t += dt;       
       printf("t=%f   x=%f  y=%f\n", t, U.x, U.y);
   }
}

Здесь мы бросаем камень со скоростью 1 м/с по горизонтали и 2 м/с по вертикали. Интервал времени равен 0.05 секунды. Программа печатает время, прошедшее со старта камня и его координаты x, y.

При заданных условиях программа напишет следующее:

t=0.050000   x=0.050000  y=0.087738
t=0.100000   x=0.100000  y=0.150950
t=0.150000   x=0.150000  y=0.189637
t=0.200000   x=0.200000  y=0.203800
t=0.250000   x=0.250000  y=0.193437
t=0.300000   x=0.300000  y=0.158550
t=0.350000   x=0.350000  y=0.099137
t=0.400000   x=0.400000  y=0.015200
t=0.450000   x=0.450000  y=-0.093263

Мы видим, что за 0.4 секунды камень пролетит на 0.4 метра вперед и упадет на землю (в следующий момент времени его координата y является отрицательной, что означает. что камень упал на землю). Максимальная высота подъема камня над землей составит приблизительно 0.2 метра. Эти результаты в очень хорошо совпадают с расчетными.

Надеюсь, вам все стало понятно и вы сможете применить эти методы в своих расчетах.