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

Задача Коши для дифференциального уравнения состоит в нахождении функции y(t), удовлетворяющей дифференциальному уравнениюn - гопорядка

y⁽ⁿ⁾ = F(t,y,y',…,y^(n-1))

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

y(t₀) = u₀, y'(t₀) = u₁,…, y^(n-1)(t₀) = u_n-1 приt = t₀.

Самым простым случаем является дифференциальное уравнение 1- го порядка, разрешенное относительно производной. Оно имеет видy' = f(t,y).Если переменнаяyявляется вектором, то речь идет о системе дифференциальных уравнений. Дифференциальные уравнения и системы порядка2и выше обычно можно свести к системам уравнений1- го порядка, введя такое количество вспомогательных функций, которое соответствует порядку уравнения. Например, чтобы свести к системе уравнений1 - го порядка уравнениеy'' = F(t,y,y'), нужно использовать подстановку:y₁ = y, y₂ = y'.В результате получим систему дифференциальных уравнений

Для численного решения обыкновенных дифференциальных уравнений произвольного порядка и систем с начальными условиями, т. е. задачи Коши, в MATLAB имеютя необходимые решатели (солверы) –ode23,ode45,ode113,ode15s,ode23s,ode23t,ode23tb. Справку по ним можно получить с помощью командыdoc ode45. Солверы с суффиксомsпредназначены для решения так называемыхжестких систем. К ним относится уравнение Ван-дер-Поляy''-μ(1 - y²)y'+y = 0, решение которого с помощьюode15sприводится во встроенной вMATLABсправочной системеHelp(приложение 1, рис. П.4).

Для всех остальных систем наиболее употребительным является ode45, реализующий алгоритм Рунге-Кутта4 - 5порядков (разные порядки точности используются для контроля шага интегрирования).

Все вышеупомянутые солверы позволяют задать дополнительный параметр options, контролирующий вычислительный процесс. Пример контроля точности вычислений рассмотрен в разделе 6.4.

В самом простом случае синтаксис вызова перечисленных выше солверов

[t,Y]=solver('fun', [t0, tn], Y0).

Здесь приняты следующие обозначения:

solver– название соответствующего солвера;

'fun'– имя файл-функции, которая вычисляет вектор-столбец правых частей системы дифференциальных уравнений;

t– вектор-столбец, содержащий значения независимой переменной (обычно это значения времени);

Y– матрица значений неизвестных функций в соответствующие моменты времени;

[t0, tn]– вектор с начальным и конечным временем наблюдения;

Y0– скаляр или вектор-столбец, в котором задаются начальные условия.

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

приведение дифференциального уравнения к системе дифференциальных уравнений первого порядка;

оформление в виде файл - функции правой части системы уравнений;

вызов подходящего солвера;

визуализация результата.

Пример:

Решить задачу Коши для дифференциального уравнения 2-го порядка

y''+4y'+3y = cost, y(0) = 1, y'(0) = 0.

Решение:

Выполнив подстановку y₁ = y, y₂ = y',получим систему дифференциальных уравнений

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

functionF=osc(t,y)

F=[y(2);-4*y(2)-3*y(1)+cos(t)];

Для вычисления решения системы на интервале [0;15]используем солверode45. С учетом начальных условий обращение к функцииode45будет иметь такой вид:

>> [T,Y]=ode45('osc',[0,15],[1;0]);

В силу проделанных замен y₁ = y, y₂ = y',первый столбец матрицыY содержит как раз значения неизвестной функцииy(t),входящей в исходное дифференциальное уравнение, а второй столбец – значения ее производной.

Отобразим на рис. 6.5 график решения и график производной от этого решения (т. е. графики изменения координаты точки и ее скорости в зависимости от времени) с помощью команды

>> plot(T,Y(:,1),T,Y(:,2),'k--')

В разделе 7.14 найдено точное (аналитическое) решение

y = cost+sint - e ^-3t+e ^-t

данного дифференциального уравнения. Выведем на рис. 6.5 график точного решения в с маркерами виде кружков:

>> t=0:15;hold

>> y=cos(t)/10+sin(t)/5-7*exp(-3*t)/20+5*exp(-t)/4;

>> plot(t,y,'ko'), legend('y(t)','dy(t)/dt','analitic solution'), grid

Рис. 6.5

Как видно из рис. 6.5, график приближенного решения совпадает с графиком точного решения.

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

Об этих солверах можно узнать во встроенной в пакет MATLABсправочной системеHelp(приложение 1, рис. П.5).