Научитесь кодировать алгоритм Рунге-Кутты 4-го порядка для решения обыкновенных дифференциальных уравнений.

Обыкновенные дифференциальные уравнения (ОДУ) распространены во всех аспектах техники и науки. Они используются для моделирования физических явлений, таких как турбулентный поток воздуха, орбиты спутников, теплопередача, вибрации и многих других. Они лежат в основе того, как люди понимают физику и математику. Хотя мы можем составить эти уравнения на основе интуиции и нашего фундаментального понимания того, как устроен мир, решить их может быть еще сложнее.

Обыкновенные дифференциальные уравнения

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

Здесь вектор y может быть пространственным состоянием, массой, температурой и т. д., а t — временем. . К сожалению, большинство ОДУ невозможно решить аналитически; это означает, что функция не может быть получена, чтобы обеспечить точное решение дифференциального уравнения. Другой способ решения этих дифференциальных уравнений — использование методов численного интегрирования.

Численная интеграция

Термин числовая интеграция был впервые использован в 1915 году, но его преимущества не были по-настоящему видны до появления современных компьютеров. Численное интегрирование — это метод аппроксимации изменения функции y во времени путем знания дифференциальных уравнений, определяющих изменение y. вовремя. Они являются оценочными, как указано, поэтому они настолько хороши, насколько хороши используемые методы и приемы. Одной из самых популярных форм численного интегрирования является метод Рунге-Кутты 4-го порядка (или RK4).

RK4 можно описать уравнениями и диаграммой ниже. Чаще всего вам придется решать ОДУ в векторной форме, поэтому показана векторная форма RK4. В уравнениях значения k представляют собой оценки наклона y, рассчитанные с использованием дифференциальных уравнений в точках, показанных на диаграмма. После определения значений k они усредняются и умножаются на временной шаг ч. Затем это значение добавляется к текущему значению y, чтобы получить следующее приближенное значение временного шага y. Этот процесс повторяется до тех пор, пока не будет достигнуто желаемое количество временных шагов. Обратите внимание, что чем меньше используемое значение h, тем точнее результат будет соответствовать истинному значению.

Скрипт MATLAB RK4

Для этого примера предположим, что у нас есть частица, движение которой описывается следующими обыкновенными дифференциальными уравнениями:

Этот набор ОДУ удобен тем, что его можно решить аналитически (решение показано ниже). Они были выбраны, чтобы наш алгоритм RK4 можно было сравнить с реальным решением. Имейте в виду, что при использовании ваших собственных ОДУ они могут быть неразрешимы аналитически, поэтому вам может не с чем сравнивать ваши результаты.

Мы можем начать код, определив функцию model, которая будет возвращать значение ОДУ, когда значения t, v и w. Убедитесь, что вы поместили это определение функции в конец вашего скрипта.

function dydt = model(t,y)
    v = y(1);
    w = y(2);
    dv_dt = w;
    dw_dt = 6*v-w;
    dydt = [dv_dt,dw_dt];
end

Теперь запустим основной скрипт. Я всегда начинаю его с очистки командного окна, очистки переменных и закрытия фигур.

clc
clear variables
close all

Чтобы использовать метод RK4, нам нужно определить начальные условия, y0, и временной массив, t. Здесь t определяется как изменяющийся от 0 до 1 с шагом 1000, а начальное условие для v равно 3. а w равно 1. Массив времени используется для создания h и позже, точные значения для v и w.

% Initial Conditions and Simulation Time
y0 = [3, 1];  % y0 = [v0, w0]
t = linspace(0,1,1000)';

Затем мы можем создать массив yₖᵤₜₜₐ для хранения приближений RK4 на каждом временном шаге. Используя начальные условия для v и w, мы можем определить первый индекс ( i = 1, t =0) массива RK4. Мы также можем определить переменную временного шага, h. Затем мы вводим цикл for для перебора длины массива времени. На каждой итерации мы будем добавлять новое значение к yₖᵤₜₜₐ, а затем использовать эту итерацию для создания другого значения для yₖᵤₜₜₐ. . Этот процесс будет повторяться, пока мы не достигнем последнего временного шага (t = 1).

% Runge-Kutta 4th-Order Algorithm
y_Kutta = zeros(length(t), 2);
y_Kutta(1, :) = y0;
h = t(2)-t(1);  % Constant time step
for i = 2:length(t)
    k1 = model(t(i-1), y_Kutta(i-1, :));
    k2 = model(t(i-1)+h/2, y_Kutta(i-1, :)+k1*h/2);
    k3 = model(t(i-1)+h/2, y_Kutta(i-1, :)+k2*h/2);
    k4 = model(t(i-1)+h, y_Kutta(i-1, :)+k3*h);
    y_Kutta(i, :) = y_Kutta(i-1, :)+(k1/6+k2/3+k3/3+k4/6)*h;
end

Как упоминалось ранее, этот набор ОДУ имеет аналитическое решение, поэтому давайте создадим массив этих значений и сравним его с нашим результатом из RK4. Для сравнения мы можем просто вычислить разницу между точным решением и нашим приближенным решением.

% Exact
y_Exact = 2*exp(2*t).*[1, 2]-exp(-3.*t)*[-1, 3];
% Calculating the Difference Between Exact and RK4 Solutions
diff_Kutta = y_Exact-y_Kutta;

Наконец, у нас есть все данные, необходимые для сравнения метода Рунге-Кутты 4-го порядка с точным решением. Последним шагом является построение результатов.

% Comparison Subplot
figure; subplot(2,1,1);  hold on;
plot(y_Exact(:, 1), y_Exact(:, 2))
plot(y_Kutta(:, 1), y_Kutta(:, 2))
title('RK4 Compared to the Exact Solution', 'Interpreter', 'Latex')
xlabel('v', 'Interpreter', 'Latex')
ylabel('w', 'Interpreter', 'Latex')
legend('Exact', 'RK4', 'Location', 'Best')
grid on; hold off
% Difference Subplot
subplot(2,1,2); hold on;
plot(t, diff_Kutta(:, 1))
plot(t, diff_Kutta(:, 2))
xlabel('t', 'Interpreter', 'Latex')
ylabel('Difference', 'Interpreter', 'Latex')
legend('RK4 [v]', 'RK4 [w]', 'Location', 'Best')
grid on; hold off

Результаты:

Как видите, результаты метода Рунге-Кутты 4-го порядка для этого примера очень точны. Результаты RK4 неотличимы от точного решения на верхнем графике, а разница между результатами для v и w и точное решение практически ничтожно. Рунге-Кутта может быть расширен за пределы 4-го порядка, чтобы уменьшить ошибку, но это, вероятно, не обязательно для большинства приложений.

Спасибо за прочтение статьи! Надеюсь, вы немного узнали о RK4! Подпишитесь на меня и ознакомьтесь с другими моими статьями о Python, математике и орбитальной механике! Если у вас есть какие-либо комментарии или проблемы, дайте мне знать!