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

Задача Коши для дифференциального уравнения состоит в нахождении функции 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, по которой производится дифференцирование и вектор, размер которого равен числу неизвестных функций системы) и один выходной аргумент (вектор правой части системы). Для нашего примера текст файл-функции будет таким:

function F=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).