Численное решение обыкновенных дифференциальных уравнений в системе MATLAB

Выработано много удобных и относительно экономных способов проводить подоб­ные вычисления. Их подробно изучают в курсе "Вычислительные методы". В системе MATLAB существует целый пакет процедур, предназначен­ных для решения систем обыкновенных дифференциальных уравнений, а если фор­мулировать еще более точно, то для решения задачи Коши для системы обыкно­венных дифференциальных уравнений (ОДУ). Причем этот пакет обеспечивает на выбор разные алгоритмы решения такой задачи, но для экономии уси­лий и облегчения работы обращение и написание требуемого дополнительного кода одинаково и не зависит от алгоритма. Поэтому рассмотрим здесь метод подготов­ки М-файлов для решения таких систем независимо от используемого алгоритма решения.

Из приведеных выше примеров понятно, что большой класс ОДУ, а именно уравнения с одной независимой переменной, которую чаще всего именуют време­нем t, разрешенные относительно старшей производной могут быть сведены к си­стеме дифференциальных уравнений первого порядка вида

y' (t) =F (t,y (t))

с начальными условиями

y(t₀) = y₀-

Если правые части соответствующей системы достаточно гладки, то описанная система имеет единственное решение, которое может быть получено численно с помощью одного из алгоритмов, используемых в системе MATLAB.

Подготовка М-файла для решения ОДУ (пример)

Рассмотрим пример подготовки М-файла для решения уравнения Ван-дер-Поля

Сделав очевидные замены у ₁ = у, y ₂ = у ', можно переписать это уравнение в виде системы двух уравнений:

Для решения этой системы уравнений необходимо написать М-файл, который опи­сывал бы правую часть системы уравнений. Для нашего случая эта функция (у ₁ и y ₂ становятся элементами двумерного вектора у(1) и у(2)) должна вы­глядеть следующим образом:

function dy = vdpl(t,y,mu)

dy = [y(2);

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

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

Для решения системы уравнений на временном интервале [0 20] и с начальными значениями y1(0) = 2, y2(0) = 0 необходимо написать

[T, Y]= ode45(@vdpl, [0 20],[2;0],[],1);

В результате решения получим вектор-столбец Т, содержащий точки на заданном временном интервале, и матрицу Y, каждая строка которой соответствует времени из Т. Для того чтобы увидеть результат решения, можно написать такую последо­вательность команд

% Вывод результатов решения

plot(T,Y(:,1), '-', T,Y(:,2), '--');

title('Решение ур-ния Ван-дер-Поля для mu =1')

xlabel ('Время Т');

ylabel('Решение Y');

legend ('y', 'y1');

Методы решения ОДУ

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

ode45 базируется на явном методе Рунге-Кутта. Это одношаговый алгоритм для вычисления y (t _n) необходимо знание решения в одной предыдущей точке y (t_n-1). Эта функция наиболее удобна для первого, 'пристрелочного' реше­ния большинства задач.

ode23 тоже базируется на явном методе Рунге-Кутта. но меньшего порядка, по­этому бывает более подходящим для получения более грубого решения (с меньшей точностью) и при наличии небольшой жесткости. Также является одношаговым методом.

ode113 использует метод переменного порядка Адамса-Бэшфорта-Милтона. Он может оказаться более эффективным, чем метод ode45, особенно при высо­ких точностях и при сложности вычисления правых частей уравнений. Метод многошаговый, поэтому для начала решения необходимо знание решения в нескольких начальных точках.

Для решения жестких систем уравнений в системе MATLAB предусмотрены 4 функции: odel5s, ode23s, ode23t, ode23tb.