Обыкновенные дифференциальные уравнения

В системе MatLab предусмотрены специальные средства решения задачи Коши для системобыкновенных дифференциальных уравнений, заданных как в явной форме = F(t,x), так и в неявной M=F(t,x), где М- матрица, - т.н. решатель ОДУ (solver ODE), обеспечивающий пользователю возможность выбора метода, задания начальных условий и др.

В простейшем варианте достаточно воспользоваться командой [T,X]=solver('F', [DT], X0,...), где DT - диапазон интегрирования, X0 - вектор начальных значений, F - имя функции вычисления правых частей системы, solver - имя используемой функции (ode45 - метод Рунге-Кутта 4 и 5-го порядков, ode23 - тот же метод 2 и 3-го поряд-ков, ode113 - метод Адамса - для нежестких систем, ode23s, ode15s - для жестких систем и др.). Версии решателя различаются используемыми методами (по умолчанию относительная погрешность 10^-3 и абсолютная 10^-6) и соответственно временем и успешностью решения. Под жесткостью здесь понимается повышенное требование к точности - использование минимального шага во всей области интегрирования. При отсутствии информации о жесткости рекомендуется попытаться получить решение посредством ode45 и затем ode15s. Если диапазон DT задан начальным и конечным значением [t_o, t_k], то количество элементов в массиве Т (и в массиве решений Х) определяется необходимым для обеспечения точности шагом; при задании DT в виде [t_o, t₁, t₂,..., t_k] или [t_o:Dt: t_k] - указанными значениями.

Например, в простейшем варианте решение уравнения =tЧe^-t в интервале t О[0, 0.5] с начальным условием x(t=0)=1:

Для иллюстрации решения системы и ряда нестандартных возможностей рассмотрим задачу выравнивания цен по уровню актива в следующей постановке.

Предположим, что изменение уровня актива у пропорционально разности между предложением s и спросом р, т.е. y'=k(s-d), k>0, и что изменение цены z пропорционально отклонению актива у от некоторого уровня y₀, т.е. z'=-m(y-y0), m>0. Естественно, что предложение и спрос зависят от цены, например, s(z)=az+s0, d(z)= d0 - cz. Соответственно возникает система дифференциальных уравнений

y' = kЧ(s(z)-d(z))
z' = - mЧ (y-y0).

Вычисление правых частей оформляем функцией:

Если выполнить решение при y₀ = 19, z₀ =2

будет выведена таблица значений искомых функций:

и их графики (рис. 8.6).

Эта система уравнений относится к числу т.н. автономных (или динамических), ибо независимая переменная в нее явно не входит; соответственно может быть установлена связь между найденными решениями: в параметрическом задании линия y=y(t), z=z(t) определяет фазовую кривую (траекторию) системы - гладкую кривую без самопересечений, замкнутую кривую или точку, которая позволяет судить об устойчивости системы.

Так, установив опции к построению двумерного фазового портрета (функция odephas2)и номера переменных состояния:
>> opt=odeset('OutputSel',[1 2], 'OutputFcn','odephas2');
>> [T,Y]=ode45('odu2', [0:0.3:9],[19 2],opt);
получаем фазовый портрет системы, свидетельствующий о ее устойчивости - гармонии между активом и ценами (рис.8.7).

Рис.8.6

Рис.8.7

В качестве другого примера подобных задач рассмотрим известную задачу динамики популяций, где рассматривается модель взаимодействия "жертв" и "хищников", в которой учитывается уменьшение численности представителей одной стороны с ростом численности другой. Модель была создана для биологических систем, но с определенными корректурами применима к конкуренции фирм, строительству финансовых пирамид, росту народонаселения, экологической проблематике и др.

Эта модель Вольтерра-Лотка с логистической поправкой описывается системой уравнений

с условиями заданной численности "жертв" и "хищников" в начальный момент t=0.

Решая эту задачу при различных значениях a, получаем различные фазовые портреты (обычный колебательный процесс и постепенная гибель популяций).

>>opt=odeset('OutputSel',[1 2], 'OutputFcn', 'odephas2');
>> [T,X]=ode45('VolterraLog', [0 10],[3 1],opt);

Рис.8.8 (a =0)

Рис.8.9 (a =0.1)

Имеется возможность построения и трехмерного фазового портрета с помощью функции odephas3. Например, решение задачи Эйлера свободного движения твердого тела:

=x₂x₃,

=-x₁x₃,

=-0.51Чx1x₂

выступает в виде:

Рис.8.10.

Рис.8.11.

наверх

следующая глава