ode15s

Решение жестких дифференциальных уравнений и ДАУ — метод переменного порядка точности

Синтаксис

[t,y] = ode15s(odefun,tspan,y0)
[t,y] = ode15s(odefun,tspan,y0,options)
[t,y,te,ye,ie] = ode15s(odefun,tspan,y0,options)
sol = ode15s(___)

Описание

пример

[t,y] = ode15s(odefun,tspan,y0), где tspan = [t0 tf], интегрирует систему дифференциальных уравнений y'=f(t,y) от t0 до tf с начальными условиями y0. Каждая строка в массиве решения y соответствует значению, возвращенному в вектор-столбце t.

Все решатели MATLAB® ODE могут решить системы уравнений формы y'=f(t,y), или проблемы, которые включают большую матрицу, M(t,y)y'=f(t,y). Решатели все использование подобные синтаксисы. Решатель ode23s только может решить проблемы с большой матрицей, если большая матрица является постоянной. ode15s и ode23t могут решить проблемы с большой матрицей, которая сингулярна, известна как дифференциально-алгебраические уравнения (ДАУ). Задайте большую матрицу с помощью опции Mass odeset.

пример

[t,y] = ode15s(odefun,tspan,y0,options) также использует настройки интегрирования, заданные options, который является созданным использованием аргумента функции odeset. Например, используйте AbsTol и опции RelTol, чтобы задать допуски абсолютной и относительной погрешности или опцию Mass, чтобы обеспечить большую матрицу.

[t,y,te,ye,ie] = ode15s(odefun,tspan,y0,options) дополнительно находит, где функции (t, y), вызвал функции события, являются нулем. В выводе te является временем события, ye является решением во время события, и ie является индексом инициированного события.

Для каждой функции события задайте, должно ли интегрирование остановиться в нуле и имеет ли направление нулевого пересечения значение. Сделайте это путем установки свойства 'Events' на функцию, такую как myEventFcn или @myEventFcn, и создания соответствующей функции: [value, isterminal, direction] = myEventFcn (t, y). Для получения дополнительной информации смотрите Местоположение События ОДУ.

пример

sol = ode15s(___) возвращает структуру, которую можно использовать с deval, чтобы оценить решение в любой точке на интервале [t0 tf]. Можно использовать любую из комбинаций входных аргументов в предыдущих синтаксисах.

Примеры

свернуть все

Простые ОДУ, которые имеют компонент единого решения, могут быть заданы как анонимная функция в вызове решателя. Анонимная функция должна принять два входных параметров (t,y), даже если одни из входных параметров не используются.

Решите ОДУ

y=-10t.

Используйте временной интервал [0,2] и начального условия y0 = 1.

tspan = [0 2];
y0 = 1;
[t,y] = ode15s(@(t,y) -10*t, tspan, y0);

Постройте решение.

plot(t,y,'-o')

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

Система уравнений:

Начальные условия и. Функциональный vdp1000 поставляется с MATLAB® и кодирует уравнения.

function dydt = vdp1000(t,y)
%VDP1000  Evaluate the van der Pol ODEs for mu = 1000.
%
%   See also ODE15S, ODE23S, ODE23T, ODE23TB.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2014 The MathWorks, Inc.

dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];

Решение этой системы с помощью ode45 с допусками относительной и абсолютной погрешности по умолчанию (1e-3 и 1e-6, соответственно) является чрезвычайно медленным, требуя, чтобы несколько минут решили и построили решение. ode45 требует, чтобы миллионы временных шагов завершили интегрирование, из-за областей жесткости, где это изо всех сил пытается соответствовать допускам.

Это - график решения, полученного ode45, который занимает много времени, чтобы вычислить. Заметьте огромное количество временных шагов, требуемых проходить через области жесткости.

Решите жесткую систему с помощью решателя ode15s, и затем постройте первый столбец решения y против моментов времени t. Решатель ode15s проходит через жесткие области с гораздо меньшим количеством шагов, чем ode45.

[t,y] = ode15s(@vdp1000,[0 3000],[2 0]);
plot(t,y(:,1),'-o')

ode15s только работает с функциями, которые используют два входных параметра, t и y. Однако можно передать в дополнительных параметрах путем определения их вне функции и передачи их в том, когда вы задаете указатель на функцию.

Решите ОДУ

При перезаписи уравнения, когда уступает система первого порядка

odefcn.m представляет эту систему уравнений как функцию, которая принимает четыре входных параметра: t, y, A и B.

function dydt = odefcn(t,y,A,B)
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = (A/B)*t.*y(1);

Решите ОДУ с помощью ode15s. Задайте указатель на функцию, таким образом, что он передает в предопределенных значениях для A и B к odefcn.

A = 1;
B = 2;
tspan = [0 5];
y0 = [0 0.01];
[t,y] = ode15s(@(t,y) odefcn(t,y,A,B), tspan, y0);

Постройте график результатов.

plot(t,y(:,1),'-o',t,y(:,2),'-.')

Решатель ode15s является хорошим предпочтительным вариантом для самых жестких проблем. Однако другие жесткие решатели могут быть более эффективными для определенных типов проблем. Этот пример решает жесткое тестовое уравнение с помощью всех четырех жестких решателей ОДУ.

Рассмотрите тестовое уравнение

y=-λy.

Уравнение становится все больше жестким как значение λ увеличения. Использование λ=1×109 и начальное условие y(0)=1 по временному интервалу [0 0.5]. С этими значениями проблема достаточно жестка, что ode45 и ode23 изо всех сил пытаются интегрировать уравнение. Кроме того, используйте odeset, чтобы передать в постоянном якобиане J=fy=-λ и включите отображение статистики решателя.

lambda = 1e9;
y0 = 1;
tspan = [0 0.5];
opts = odeset('Jacobian',-lambda,'Stats','on');

Решите уравнение с ode15s, ode23s, ode23t и ode23tb. Сделайте подграфики для сравнения.

subplot(2,2,1)
tic, ode15s(@(t,y) -lambda*y, tspan, y0, opts), toc
104 successful steps
1 failed attempts
212 function evaluations
0 partial derivatives
21 LU decompositions
210 solutions of linear systems
Elapsed time is 1.270468 seconds.
title('ode15s')
subplot(2,2,2)
tic, ode23s(@(t,y) -lambda*y, tspan, y0, opts), toc
63 successful steps
0 failed attempts
191 function evaluations
0 partial derivatives
63 LU decompositions
189 solutions of linear systems
Elapsed time is 0.254141 seconds.
title('ode23s')
subplot(2,2,3)
tic, ode23t(@(t,y) -lambda*y, tspan, y0, opts), toc
95 successful steps
0 failed attempts
125 function evaluations
0 partial derivatives
28 LU decompositions
123 solutions of linear systems
Elapsed time is 0.397964 seconds.
title('ode23t')
subplot(2,2,4)
tic, ode23tb(@(t,y) -lambda*y, tspan, y0, opts), toc
71 successful steps
0 failed attempts
167 function evaluations
0 partial derivatives
23 LU decompositions
236 solutions of linear systems
Elapsed time is 0.314648 seconds.
title('ode23tb')

Жесткие решатели, которые все выполняют хорошо, но ode23s завершает интеграцию с наименьшим количеством шагов и запускает самое быстрое для этой конкретной проблемы. Поскольку постоянный якобиан задан, ни один из решателей не должен вычислять частные производные, чтобы вычислить решение. Определение якобиана приносит пользу ode23s больше всего, поскольку это обычно оценивает якобиан на каждом шаге.

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

Уравнение Ван дер Поля является ОДУ второго порядка

y1-μ(1-y12)y1+y1=0.

Решите уравнение Ван дер Поля с μ=1000 использование ode15s. Функциональный vdp1000.m поставляется с MATLAB® и кодирует уравнения. Задайте один вывод, чтобы возвратить структуру, содержащую информацию о решении, таком как точки оценки и решатель.

tspan = [0 3000];
y0 = [2 0];
sol = ode15s(@vdp1000,tspan,y0)
sol = struct with fields:
     solver: 'ode15s'
    extdata: [1x1 struct]
          x: [1x592 double]
          y: [2x592 double]
      stats: [1x1 struct]
      idata: [1x1 struct]

Используйте linspace, чтобы сгенерировать 2 500 точек в интервале [0 3000]. Оцените первый компонент решения в этих точках с помощью deval.

x = linspace(0,3000,2500);
y = deval(sol,x,1);

Постройте решение.

plot(x,y)

Расширьте решение tf=4000 использование odextend и добавляет результат в исходный график.

tf = 4000;
sol_new = odextend(sol,@vdp1000,tf);
x = linspace(3000,tf,350);
y = deval(sol_new,x,1);
hold on
plot(x,y,'r')

Этот пример повторно формулирует систему ОДУ как система дифференциальных алгебраических уравнений (ДАУ). Задачей Робертсона, найденной в hb1ode.m, является классическая тестовая задача для программ, которые решают жесткие ОДУ. Система уравнений

hb1ode решает эту систему ОДУ к устойчивому состоянию с начальными условиями, и. Но уравнения также удовлетворяют линейный закон сохранения,

С точки зрения решения и начальных условий, закон сохранения

Система уравнений может быть переписана как система ДАУ при помощи закона сохранения, чтобы определить состояние. Это повторно формулирует проблему как систему ДАУ

Дифференциальный индекс этой системы равняется 1, поскольку только одна производная требуется, чтобы делать это системой ОДУ. Поэтому никакие дальнейшие преобразования не требуются прежде, чем решить систему.

Функциональный robertsdae кодирует эту систему ДАУ. Сохраните robertsdae.m в своей текущей папке, чтобы запустить пример.

function out = robertsdae(t,y)
out = [-0.04*y(1) + 1e4*y(2).*y(3)
   0.04*y(1) - 1e4*y(2).*y(3) - 3e7*y(2).^2
   y(1) + y(2) + y(3) - 1 ];

Полный пример кода для этой формулировки задачи Робертсона доступен в hb1dae.m.

Решите систему ДАУ с помощью ode15s. Сопоставимые начальные условия для y0 очевидны на основе закона сохранения. Используйте odeset, чтобы установить опции:

  • Используйте постоянную большую матрицу, чтобы представлять левую сторону системы уравнений.

  • Установите допуск относительной погрешности на 1e-4.

  • Используйте абсолютный допуск 1e-10 для второго компонента решения, поскольку шкала отличается существенно от других компонентов.

  • Оставьте параметр «MassSingular» по умолчанию «maybe», чтобы проверить автоматическое обнаружение ДАУ.

y0 = [1; 0; 0];
tspan = [0 4*logspace(-6,6)];
M = [1 0 0; 0 1 0; 0 0 0];
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6]);
[t,y] = ode15s(@robertsdae,tspan,y0,options);

Постройте решение.

y(:,2) = 1e4*y(:,2);
semilogx(t,y);
ylabel('1e4 * y(:,2)');
title('Robertson DAE problem with a Conservation Law, solved by ODE15S');

Входные параметры

свернуть все

Функции, чтобы решить, заданный как указатель на функцию, который задает функции, которые будут интегрированы.

Функциональный dydt = odefun(t,y), для скалярного t и вектор-столбца y, должен возвратить вектор-столбец dydt типа данных single или double, который соответствует f(t,y). odefun должен принять оба входных параметра, t и y, даже если один из аргументов не используется в функции.

Например, чтобы решить y'=5y3, используйте функцию:

function dydt = odefun(t,y)
dydt = 5*y-3;

Для системы уравнений вывод odefun является вектором. Каждый элемент в векторе является решением одного уравнения. Например, чтобы решить

y'1=y1+2y2y'2=3y1+2y2

используйте функцию:

function dydt = odefun(t,y)
dydt = zeros(2,1);
dydt(1) = y(1)+2*y(2);
dydt(2) = 3*y(1)+2*y(2);

Для получения информации о том, как предоставить дополнительные параметры функциональному odefun, смотрите Функции Параметризации.

Пример: @myFcn

Типы данных: function_handle

Интервал интегрирования, заданного как вектор. В минимуме tspan должен быть двумя векторами элемента [t0 tf], задающий начальные и итоговые времена. Чтобы получить решения в конкретные моменты времени между t0 и tf, используйте более длинный вектор формы [t0,t1,t2,...,tf]. Элементы в tspan должны все увеличиваться или все уменьшение.

Решатель налагает начальные условия, данные y0 в начальное время tspan(1), затем объединяется от tspan(1) до tspan(end):

  • Если tspan имеет два элемента, [t0 tf], то решатель возвращает решение, оцененное в каждом внутреннем этапе интеграции в интервале.

  • Если tspan имеет больше чем два элемента [t0,t1,t2,...,tf], то решатель возвращает решение, оцененное в данных точках. Однако решатель не продвигается точно в каждую точку, заданную в tspan. Вместо этого решатель использует свои собственные внутренние шаги, чтобы вычислить решение, затем оценивает решение в требуемых точках в tspan. Решения, произведенные в заданных точках, имеют тот же порядок точности как решения, вычисленные на каждом внутреннем шаге.

    Определение нескольких промежуточных точек имеет мало эффекта на эффективность вычисления, но для больших систем это может влиять на управление памятью.

Значения tspan используются решателем, чтобы вычислить подходящие значения для InitialStep и MaxStep:

  • Если tspan содержит несколько промежуточных точек [t0,t1,t2,...,tf], то заданные точки дают индикацию относительно шкалы для проблемы, которая может влиять на значение InitialStep, используемого решателем. Поэтому решение, полученное решателем, может отличаться в зависимости от того, задаете ли вы tspan как двухэлементный вектор или как вектор с промежуточными точками.

  • Начальные и окончательные значения в tspan используются, чтобы вычислить максимальный размер шага MaxStep. Поэтому изменение начальных или окончательных значений в tspan могло привести к решателю с помощью различной последовательности шага, которая может изменить решение.

Пример: [1 10]

Пример: [1 3 5 7 9 10]

Типы данных: single | double

Начальные условия, заданные как вектор. y0 должен быть той же длиной как векторный вывод odefun, так, чтобы y0 содержал начальное условие для каждого уравнения, определенного в odefun.

Типы данных: single | double

Структура опции, заданная как массив структур. Используйте функцию odeset, чтобы создать или изменить структуру опций. См. Сводные данные Опций ОДУ для списка опций, совместимых с каждым решателем.

Пример: options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot) задает допуск относительной погрешности 1e-5, включает отображение статистики решателя и задает выходную функцию @odeplot, чтобы построить решение, когда это вычисляется.

Типы данных: struct

Выходные аргументы

свернуть все

Точки оценки, возвращенные как вектор-столбец.

  • Если tspan содержит два элемента, [t0 tf], то t содержит внутренние точки оценки, раньше выполнял интегрирование.

  • Если tspan содержит больше чем два элемента, то t совпадает с tspan.

Решения, возвращенные как массив. Каждая строка в y соответствует решению в значении, возвращенном в соответствующей строке t.

Время событий, возвращенных как вектор-столбец. Время наступления события в te соответствует решениям, возвращаемым в ye, и ie определяет, какое событие произошло.

Решение во время событий, возвращенных как массив. Время наступления события в te соответствует решениям, возвращаемым в ye, и ie определяет, какое событие произошло.

Индекс исчезающей функции события, возвращенной как вектор-столбец. Время наступления события в te соответствует решениям, возвращаемым в ye, и ie определяет, какое событие произошло.

Структура для оценки, возвращенной как массив структур. Используйте эту структуру с функцией deval, чтобы оценить решение в любой точке в интервале [t0 tf]. Массив структур sol всегда включает эти поля:

Поле структурыОписание

sol.x

Вектор строка с шагом, выбранным решателем.

sol.y

Решения. Каждый столбец sol.y(:,i) содержит решение во время sol.x(i).

sol.solver

Имя решателя.

Кроме того, если вы задаете опцию Events, и события обнаруживаются, затем sol также включает эти поля:

Поле структурыОписание

sol.xe

Точки, когда события имели место. sol.xe(end) содержит точное место терминального события, если таковые имеются.

sol.ye

Решения, которые соответствуют событиям в sol.xe.

sol.ie

Индексы в вектор, возвращенный функцией, заданы в опции Events. Значения указывают, какое событие решатель обнаружил.

Алгоритмы

ode15s является переменным шагом, переменный порядок (VSVO) решатель на основе числовых формул дифференцирования (NDFs) порядков 1 - 5. Опционально, это может использовать формулы дифференцирования назад (BDF, также известные как метод Гира), которые обычно менее эффективны. Как ode113, ode15s является многоступенчатым решателем. Используйте ode15s, если ode45 перестал работать или очень неэффективен, и вы подозреваете, что проблема жестка, или при решении дифференциально-алгебраического уравнения (DAE) [1], [2].

Ссылки

[1] Шемпин, L. F. и М. В. Рейчелт, “Пакет ODE MATLAB”, SIAM Journal на Научных вычислениях, Издании 18, 1997, стр 1–22.

[2] Шемпин, L. F. М. В. Рейчелт и Дж.А. Кирженка, “Решая Индексные 1 ДАУ в MATLAB и Simulink”, Анализ SIAM, Издание 41, 1999, стр 538–552.

Представлено до R2006a

Для просмотра документации необходимо авторизоваться на сайте