Использование решателей систем ОДУ
В описанных далее функциях для решения систем дифференциальных уравнений приняты следующие обозначения и правила:
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}].
Решатели используют в списке различные параметры. В приведенной ниже таблице они даны для решателей обычных (в том числе жестких) дифференциальных уравнений.
|
|
|
||||||
|
|
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) помечающие их надписи.