Пошаговый метод численного интегрирования дифференциальных уравнений с использованием MATLAB
Дифференциальные уравнения описывают множество фундаментальных физических законов. Они используются для моделирования траекторий ракет, для оценки уровня заражения смертельной болезнью, для моделирования экономического роста страны и т. д. Они окружают нас повсюду и при правильном использовании могут дать впечатляющие результаты. Большинство дифференциальных уравнений, имеющих практическое применение, представляют собой обыкновенные дифференциальные уравнения (ОДУ). Это уравнения, состоящие из одной или нескольких функций (и их производных) с одной независимой переменной. Независимой переменной в большинстве ОДУ является время или t. К сожалению, большинство ОДУ слишком сложно или невозможно решить в явном виде. Следующим лучшим шагом в получении результата для сложного ОДУ является численное интегрирование.
Численное интегрирование используется для аппроксимации эволюции переменных ОДУ с течением времени. В этой статье мы сосредоточимся на проблемах начального значения. Для этих типов задач численное интегрирование использует начальное время и начальные условия, проходит через заранее определенный интервал времени и вычисляет решение для переменных на каждом временном шаге. Преимущество численного интегрирования заключается в том, что ОДУ можно аппроксимировать с высокой степенью точности в зависимости от используемого метода. Кроме того, благодаря большому количеству вычислительной мощности, доступной сегодня, эти методы могут выполняться за доли секунды.
MATLAB имеет множество числовых интеграторов и является выбором многих инженеров и физиков за его вычислительные возможности и относительно простой синтаксис. Поскольку MATLAB разработан и поддерживается MathWorks, они предоставляют подробные объяснения и списки каждого из своих решателей ОДУ. Выбор правильного решателя для вашей проблемы может иметь решающее значение, но для большинства приложений лучше всего начать с решателя ode45. Решатель ode45 — это решатель Runge Kutta 4-го и 5-го порядка. Я лучше всего учусь на примере, так что давайте перейдем к коду.
Во-первых, давайте определим, какие обыкновенные дифференциальные уравнения мы хотим изучать. Одним из самых интересных наборов ОДУ является система Лоренца. Этот набор дифференциальных уравнений был изучен Эдвардом Лоренцем и может быть использован для создания упрощенной модели атмосферной конвекции. Поскольку мы больше сосредоточимся на том, как численно решать ОДУ, если вам нужна дополнительная информация о самой модели, вы можете посетить ссылку на Википедию выше. Обыкновенные дифференциальные уравнения, описывающие систему Лоренца, приведены ниже:
Здесь нас интересуют переменные x, y и z. Мы определим константы σ, ρ и β позже в коде. Чтобы численно интегрировать эти уравнения, нам нужно определить все, что нужно функции ode45: начальные условия, временной интервал и модель. Наши начальные условия имеют форму Y ниже (наши результаты также будут иметь эту форму), а модель представляет собой производную по времени вектора Y (или ОДУ системы Лоренца). :
Начнем код с определения модели (или производной по времени от Y, dYdt) как функции MATLAB. Следующий код показывает именно это. Начнем с извлечения x, y и z из вектора Y (не путать с переменной y), который передается в model. Затем мы определяем вектор производной по времени, dYdt, который будет возвращен этой функцией. Это ОДУ, которые были определены системой Лоренца. Здесь мы определяем σ, ρ и β как 10, 20 и 1,5 соответственно.
% ODEs or Model
function dYdt = model(t, Y)
x = Y(1);
y = Y(2);
z = Y(3);
dYdt = zeros(size(Y));
dYdt(1) = 10*(y-x); % dxdt
dYdt(2) = x*(20-z)-y; % dydt
dYdt(3) = x*y-1.5*z; % dzdt
end
Функции MATLAB должны быть определены в конце кода, поэтому следующий код необходимо поместить в начало скрипта, иначе MATLAB будет недоволен вами. Во-первых, согласно стандартной практике, мы очищаем командное окно, очищаем переменные и закрываем все окна фигур. Это позволяет нам начать с чистого листа, прежде чем запускать остальную часть скрипта. Далее мы определим два других важных входа в ode45, начальные условия и временной интервал интегрирования. Они оба были выбраны произвольно для этого примера, но у вас могут быть определенные начальные условия и временной интервал, которые вас интересуют для вашей проблемы.
clc clear variables close all % Initial Conditions and Time Span Y0 = [1; 0; 0]; % [x; y; z] tspan = [0 50]; % [Time_Start Time_Final]
Следующий шаг — просто передать нашу функцию MATLAB, интересующий временной интервал и начальные условия, которые мы уже определили, в ode45. Мы также можем определить опцию opts для ode45, которая обеспечит точность нашего решения в той степени, в которой вы хотите. Здесь мы используем odeset, чтобы установить очень низкую относительную устойчивость решения. Возможно, вам не понадобится такой низкий уровень для вашей конкретной проблемы, поэтому определите, что лучше для вас. После запуска ode45 с этими аргументами вы можете извлечь решение для x, y и z.
% Numerically Integrating
opts = odeset('RelTol', 1e-13);
[t, Y] = ode45(@model, tspan, Y0, opts);
x = Y(:, 1);
y = Y(:, 2);
z = Y(:, 3);
Наконец, вы можете построить график результатов для x, y и z из численного интегрирования с помощью plot3 функция (в данном случае, поскольку у нас есть три переменные).
% Displaying Results
figure; hold on
plot3(x, y, z, 'b')
xlabel('x')
ylabel('y')
zlabel('z')
view(45, 45)
grid minor
Вы получите следующую систему Лоренца из численного интегрирования, если все прошло правильно.
Это весь код, необходимый для численного интегрирования обыкновенных дифференциальных уравнений с помощью MATLAB. Мощная библиотека решателей ОДУ в MATLAB делает этот процесс относительно простым, если вы знаете, как это сделать. Если вы заинтересованы в численной интеграции в Python, я рекомендую вам ознакомиться с этой статьей:
Вы также можете научиться создавать свой собственный числовой интегратор, используя методы Рунге-Кутты в MATLAB:
Спасибо за чтение! Пожалуйста, оставьте комментарий, если вам нужна помощь в реализации этого. Я был бы более чем счастлив помочь. Если вас интересуют другие статьи по кодированию, инженерии или космосу, подпишитесь на меня!
