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

Задача Коши для дифференциального уравнения состоит в нахождении функции 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).

Вопросы для самопроверки

1. Как осуществляется ввод полиномов в MATLAB?

2. Какие команды осуществляют умножение, деление полиномов?

3. Какая команда вычисляет корни полинома?

4. Какая команда осуществляет построение полинома по заданному вектору его корней?

5. Какая команда вычисляет значение полинома по заданному значению его аргумента?

6. Каким образом можно задать первый входной аргумент команд fzero и fsolve?

7. Какая команда MATLAB находит минимум функции одной переменной?

8. Какая команда MATLAB находит минимум функции нескольких переменных?

9. Какие команды вычисляют значения определенных интегралов?

10. Какие команды (солверы) возвращают численные решения обыкновенных дифференциальных уравнений?

Не нашли, что искали? Воспользуйтесь поиском по сайту: