Использование
решателей систем ОДУ
В описанных
далее функциях для решения систем дифференциальных уравнений приняты следующие
обозначения и правила:
-
options
— аргумент, создаваемый функцией odeset (еще одна функция — odeget или bvpget
(только для bvp4c)— позволяет вывести параметры, установленные по умолчанию
или с помощью функции odeset /bvpset);
-
tspan
— вектор, определяющий интервал интегрирования [tO tfinal]. Для получения
решений в конкретные моменты времени t0, tl,..., tfinal (расположенные в
порядке уменьшения или увеличения) нужно использовать tspan = [t0 tl ...
tfinal];
-
у0 — вектор
начальных условий;
-
pi, р2,.„
— произвольные параметры, передаваемые в функцию F;
-
Т, Y —
матрица решений Y, где каждая строка соответствует времени, возвращенном
в векторе-столбце Т.
Перейдем к
описанию функций для решения систем дифференциальных уравнений:
-
[T.Y]
= solver(@F,tspan,у0) — где вместо solver подставляем имя конкретного решателя
— интегрирует систему дифференциальных уравнений вида
у'=F(t,y)
на
интервале tspan с начальными условиями у0. @F — дескриптор ODE-функции.
Каждая строка в массиве решений Y соответствует значению времени, возвращаемому
в векторе-столбце Т;
-
[T.Y]
= solver(@F, tspan, yO. options) — дает решение, подобное описанному выше,
но с параметрами, определяемыми значениями аргумента options, созданного
функцией odeset. Обычно используемые параметры включают допустимое значение
относительной погрешности RelTol (по умолчанию le-З) и вектор допустимых
значений абсолютной погрешности AbsTol (все компоненты по умолчанию равны
1е-6);
-
[T.Y] =
solver(@F,tspan,yO,options,pl,p2...) — дает решение, подобное описанному
выше, передавая дополнительные параметры pi, р2,... в m-файл F всякий раз,
когда он вызывается. Используйте options=[], если никакие параметры не задаются;
-
[T.Y.TE.YE.IE]
= solver(@F.tspan,yO,options) — в дополнение к описанному решению содержит
свойства Events, установленные в структуре options ссылкой на функции событий.
Когда эти функции событий от (t, у, равны нулю, производятся действия в
зависимости от значения трех векторов value, isterminal, di rection (их
величины можно установить в m-файлах функций событий). Для i-й функции событий
value(i) —значение функции, isterminal (i) — прекратить интеграцию при достижении
функцией нулевого значения, direction^) = 0, если все нули функции событий
нужно вычислять (по умолчанию), +1 — только те нули, где функция событий
увеличивается, -1 — только те нули, где функция событий уменьшается. Выходной
аргумент ТЕ — вектор-столбец времен, в которые происходят события (events),
строки YE являются соответствующими решениями, а индексы в векторе IE определяют,
какая из i функций событий (event) равна нулю в момент времени, определенный
ТЕ. Когда происходит вызов функции без выходных аргументов, по умолчанию
вызывается выходная функция odeplot для построения вычисленного решения.
В качестве альтернативы можно, например, установить свойство OutputFcn в
значение ' odephas2' или 'odephas3' для построения двумерных или трехмерных
фазовых плоскостей.
-
[T.X.Y]
= sim(@model,tspan.-y0.options,ut.p1,р2..„) — использует модель SIMULINK,
вызывая соответствующий решатель из нее. Пример:
[T.X.Y]
- sim(@model....).
Параметры
интегрирования (options) могут быть определены и в m-файле, и в командной строке
с помощью команды odeset. Если параметр определен в обоих местах, определение
в командной строке имеет приоритет.
Решатели используют
в списке параметров различные параметры:
-
NormControl
— управление ошибкой в зависимости от нормы вектора решения [on | {off}].
Установите 'on', чтобы norm(e) <= max(RelTol*norm(y), AbsTol). По умолчанию
все решатели используют более жесткое управление по каждой из составляющих
вектора решения;
-
RelTol
— относительный порог отбора [положительный скаляр]. По умолчанию 1е-3 (0.1%
точность) во всех решателях; оценка ошибки на каждом шаге интеграции e(i)
<= max(RelTol*abs(y(i)), AbsTol(i));
-
AbsTol
— абсолютная точность [положительный скаляр или вектор {1е-6}].Скаляр вводится
для всех составляющих вектора решения, а вектор указывает на компоненты
вектора решения. AbsTol по умолчанию 1е-6 во всех решателях;
-
Refine
- фактор уточнения вывода [положительное целое число] — умножает число точек
вывода на этот множитель. По умолчанию всегда равен 1, кроме ODE45, где
он 4. Не применяется, если tspan > 2;
-
OutputFcn
— дескриптор функция вывода [function] — имеет значение в том случае, если
решатель вызывается без явного указания функции вывода, OutputFcn по умолчанию
вызывает функцию odeplot. Эту установку можно поменять именно здесь;
-
OutputSel
— индексы отбора [вектор целых чисел] Установите компоненты, которые поступают
в OutputFcn. OutputSel по умолчанию выводит все компоненты;
-
Stats —
установите статистику стоимости вычислений [on {off}];
-
Jacobian
— функция матрицы Якоби [function constant matrix]. Установите это свойство
на дескриптор функции FJac (если FJac(t, у) возвращает dF/dy) или
на имя постоянной матрицы
dF/dy;
-
Jpattern
— график разреженности матрицы Якоби [имя разреженной матрицы]. Матрица
S с S(i,j) = 1, если составляющая i F(t, у) зависит от составляющей j величины
у, и 0 в противоположном случае;
-
Vectorized
— векторизованная ODE-функция [on | {off}]. Устанавливается на 'on', если
ODE-функция F F(t,[yl y2...]) возвращает вектор [F(t, yl) F(t, y2) ...];
-
Events
— [function] — введите дескрипторы функций событий, содержащих собственно
функцию в векторе value, и векторы isterminal и direction (см выше);
-
Mass —
матрица массы [constant matrix function]. Для задач М*у' = f(t, у) установите
имя постоянной матрицы. Для задач с переменной М введите дескриптор функции,
описывающей матрицу массы;
-
MstateDependence
— зависимость матрицы массы от у [none | {weak} | strong] — установите 'nоnе'
для уравнений M(t)*y' = F(t, у). И слабая ('weak'), и сильная ('strong')
зависимости означают M(t, у), но 'weak' приводит к неявным алгоритмам решения,
использующим аппроксимации при решении алгебраических уравнений;
-
MassSingular
— матрица массы М сингулярная [yes no | {maybe}] (да/нет/может быть);
-
MvPattern
— разреженность (dMv/dy), график разреженности (см функцию spy) — введите
имя разреженной матрицы S с S(i,j) = 1 для любого k, где (i, k) элемент
матрицы массы M(t, у) зависит от проекции] переменной у, и 0 в противном
случае;
-
Initial
Step — предлагаемый начальный размер шага, по умолчанию каждый решатель
определяет его автоматически по своему алгоритму;
-
Initial
SI ope — вектор начального уклона ур0 ур0 = F(t0,y0)/M(t0, y0);
-
MaxStep
— максимальный шаг, по умолчанию во всех решателях равен одной десятой интервала
tspan;
-
BDF (Backward
Differentiation Formulas) [on | {off}] — указывает, нужно ли использовать
формулы обратного дифференцирования (методы Gear) вместо формул численного
дифференцирования, используемых в ode 15s по умолчанию;
-
MaxOrder
- Максимальный порядок ODE15S [1 | 2 | 3 4 | {5}].
Решатели используют
в списке различные параметры. В приведенной ниже таблице они даны для решателей
обычных (в том числе жестких) дифференциальных уравнений.
|
Параметры
|
Ode45
|
Ode23
|
Ode11s
|
Ode15s
|
ode23s
|
|
|
RelTol,AbsTol
|
+
|
+
|
+
|
+
|
+
|
|
|
OutputFcaOutputSel,
Refine, Stats
|
+
|
+
|
+
|
+
|
+
|
|
|
Events
|
+
|
+
|
+
|
+
|
+
|
|
|
MaxStep,
InitlalStep
|
+
|
+
|
+
|
+
|
+
|
|
|
Jconstant,
Jacobl an,
|
|
|
|
|
|
|
|
Jpattern,
Vectorized
|
-
|
-
|
-
|
+
|
+
|
|
|
Mass
|
-
|
-
|
-
|
+
|
+
|
|
|
MassConstant
|
-
|
-
|
-
|
+
|
-
|
|
|
MaxOrder,
BOF
|
-
|
|
-
|
+
|
-
|
|
Решатель bvp4c
имеет очень небольшое число параметров, но можно вводить не только матрицу Якоби
интегрируемой функции, но и матрицу Якоби, содержащую частные производные функции
граничных условий по границам интервала и по неизвестным параметрам.
Покажем применение
решателя ОДУ на ставшем классическом примере — решении уравнения Ван-дер-Поля,
записанного в виде системы из двух дифференциальных уравнений:
y'
1=
y
2 ;
y'
2=
100*(1-y
1
)^2
*
y
2
-y
1
при
начальных условиях
y
1
,(0)
= 0;
y
2
(0)
= 1.
Перед решением
нужно записать систему дифференциальных уравнений в виде ode-функции. Для этого
в главном меню выберем File
>
New > M-File и введем
function
dydt = vdp100(t.y)
dydt
= zeros(2.1);
% a
column vector
dydt(l)
= y(2);
dydtC2)
= 100*(1 -у(1^)2)*у(2) -y(1);
Сохраним m-файл-функцию.
Тогда решение решателем ode15s и сопровождающий его график можно получить, используя
следующие команды:
»
[T,Y]=odel5s(@vdp100.[0 30].[2 0]);
» plot(T.Y)
»
hold on:gtext('yl').gtext('y2')
Последние
команды позволяют с помощью мыши нанести на графики решений y
1
=
y(1) и у
2
= y(2) помечающие их надписи.
Содержание раздела