Решение дифференциальных уравнений в среде MatLab

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

Решатели реализуют следующие методы решения систем дифференциальных уравнений, причем для решения жестких систем уравнений рекомендуется использовать только специальные решатели ode15s , ode23s, ode23t, ode23tb:

· ode45 - одношаговые явные методы Рунге-Кутта 4-го и 5-го порядка. Это классический метод, рекомендуемый для начальной пробы решения. Во многих случаях он дает хорошие результаты;

· ode23 - одношаговые явные методы Рунге-Кутта 2-го и 4-го порядка. При умеренной жесткости системы ОДУ и низких требованиях к точности этот мето;. может дать выигрыш в скорости решения;

· ode113 - многошаговый метод Адамса-Башворта-Мултона переменного порядка Это адаптивный метод, который может обеспечить высокую точность решения

· ode23tb - неявный метод Рунге-Кутта в начале решения и метод, использующий формулы обратного дифференцирования 2-го порядка в последующем

Несмотря на сравнительно низкую точность, этот метод может оказаться более эффективным, чем ode15s;

· ode15s - многошаговый метод переменного порядка (от 1 до 5, по умолчанию 5), использующий формулы численного дифференцирования. Это адаптивный метод, его стоит применять, если решатель ode45 не обеспечивает решения;

· ode23s - одношаговый метод, использующий модифицированную формулу Розенброка 2-го порядка. Может обеспечить высокую скорость вычислений при низкой точности решения жесткой системы дифференциальных уравнений;

· ode23t - метод трапеций с интерполяцией. Этот метод дает хорошие результаты при решении задач, описывающих колебательные системы с почти гармоническим выходным сигналом;

Все решатели могут решать системы уравнений явного вида у' = F(t, y). Решатели ode15s и ode23t способны найти корни дифференциально-алгебраических уравнений M(t)y' = F(t, у), где М называется матрицей массы. Решатели ode15s, ode23s, ode23t и ode23tb могут решать уравнения неявного вида M(t,y) у' = F(t, у).

· ode23tb, ode23s служат для решения жестких дифференциальных уравнений.

· ode15s - жестких дифференциальных и дифференциально-алгебраических уравнений,

· ode23t - умеренно жестких дифференциальных и дифференциально-алгебраических уравнений.

В описанных далее функциях для решения систем дифференциальных уравнений приняты следующие обозначения и правила:

· options - аргумент, создаваемый функцией odeset - позволяет вывести параметры, установленные по умолчанию;

· tspan - вектор, определяющий интервал интегрирования [t₀ t_final]. Для получения решений в конкретные моменты времени t₀, t_l,..., t_final (расположенные в порядке уменьшения или увеличения) нужно использовать tspan = [t₀ t_l ... t_final];

· у0 - вектор начальных условий;

· Т, Y - матрица решений Y, где каждая строка соответствует времени, возвращенном в векторе-столбце Т.

Этапы решения ОДУ:

1 Приведение дифференциального уравнения к системе дифференциальных уравнений первого порядка. Для этого вводится столько дополнительных функций, каков порядок уравнения.

2 Написание специальной файл-функции для системы уравнений. Файл-функция содержит два входных аргумента: переменную t, по которой производится дифференцирование, даже если она входит в уравнение неявно, и вектор, размер которого равен числу неизвестных функций системы.

3 Вызов подходящего солвера (встроенной функции). Входными аргументами солвера, в простом случае, являются имя файл-функции в апострофах, вектор с начальным и конечным значениями переменной

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

y'₁ = y₂ ;

y'₂ = Mu*(1 - y₁^2)*y₂ - y₁;

при начальных условиях

y₁(0) = 2;

y₂(0) = 0;

Перед решением нужно записать систему дифференциальных уравнений в виде ode-функции. Для этого в главном меню выберем File > New > M-File и введем

function dydt = vanderpol(t, y)

% функция правых частей дифф.уравнения

dydt = [y(2);

(1-y(1)^2)*y(2)-y(1)];

Сохраним m-файл-функцию. Тогда решение решателем ode45 и сопровождающий его график можно получить, используя следующие команды:

tspan = [0, 20]; % интервал интегрирования

y0 = [2; 0]; % начальные условия

[t,y] = ode45(@vanderpol, tspan, y0, []);

% Результат в виде графика

plot(t,y(:,1))

xlabel('t')

ylabel('solution y')

title('van der Pol Equation, \mu = 1')