Решите жесткие дифференциальные уравнения — метод трапеций + формула дифференцирования назад
[t,y] = ode23tb(odefun,tspan,y0)[t,y] = ode23tb(odefun,tspan,y0,options)[t,y,te,ye,ie] = ode23tb(odefun,tspan,y0,options)sol = ode23tb(___)[t,y] = ode23tb(odefun,tspan,y0)tspan = [t0 tf], интегрирует систему дифференциальных уравнений от t0 до tf с начальными условиями y0. Каждая строка в массиве решения y соответствует значению, возвращенному в векторе - столбце t.
Все решатели MATLAB® ODE могут решить системы уравнений формы , или проблемы, которые включают большую матрицу, . Решатели все использование подобные синтаксисы. Решатель ode23s только может решить проблемы с большой матрицей, если большая матрица является постоянной. ode15s и ode23t могут решить проблемы с большой матрицей, которая сингулярна, известна как дифференциально-алгебраические уравнения (ДАУ). Задайте большую матрицу с помощью опции Mass odeset.
[t,y] = ode23tb(odefun,tspan,y0,options) options, который является созданным использованием аргумента функции odeset. Например, используйте AbsTol и опции RelTol, чтобы задать допуски абсолютной и относительной погрешности или опцию Mass, чтобы обеспечить большую матрицу.
[t,y,te,ye,ie] = ode23tb(odefun,tspan,y0,options) te является временем события, ye является решением во время события, и ie является индексом инициированного события.
Для каждой функции события задайте, должно ли интегрирование отключить в нуле и имеет ли направление нулевого пересечения значение. Сделайте это путем установки свойства 'Events' на функцию, такую как myEventFcn или @myEventFcn, и создания соответствующей функции: [value, isterminal, direction] = myEventFcn (t, y). Для получения дополнительной информации смотрите Местоположение События ОДУ.
sol = ode23tb(___) deval, чтобы оценить решение в любой точке на интервале [t0 tf]. Можно использовать любую из комбинаций входных аргументов в предыдущих синтаксисах.
Простые ОДУ, которые имеют компонент единого решения, могут быть заданы как анонимная функция в вызове решателя. Анонимная функция должна принять два входных параметра (t,y), даже если один из входных параметров не используется.
Решите ОДУ
Используйте временной интервал [0,2] и начального условия y0 = 1.
tspan = [0 2]; y0 = 1; [t,y] = ode23tb(@(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, который занимает много времени, чтобы вычислить. Заметьте огромное количество временных шагов, требуемых проходить через области жесткости.
Решите жесткую систему с помощью решателя ode23tb, и затем постройте график первого столбца решения y против моментов времени t. Решатель ode23tb проходит через жесткие области с гораздо меньшим количеством шагов, чем ode45.
[t,y] = ode23tb(@vdp1000,[0 3000],[2 0]);
plot(t,y(:,1),'-o')
ode23s только работает с функциями, которые используют два входных параметра, 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);
Решите ОДУ с помощью ode23tb. Задайте указатель на функцию, таким образом, что он передает в предопределенных значениях для A и B к odefcn.
A = 1; B = 2; tspan = [0 5]; y0 = [0 0.01]; [t,y] = ode23tb(@(t,y) odefcn(t,y,A,B), tspan, y0);
Постройте график результатов.
plot(t,y(:,1),'-o',t,y(:,2),'-.')
Решатель ode15s является хорошим предпочтительным вариантом для самых жестких проблем. Однако другие жесткие решатели могут быть более эффективными для определенных типов проблем. Этот пример решает жесткое тестовое уравнение с помощью всех четырех жестких решателей ОДУ.
Рассмотрите тестовое уравнение
Уравнение становится все более и более жестким как значение
увеличений. Используйте
и начальное условие
по временному интервалу [0 0.5]. С этими значениями проблема достаточно жестка, что ode45 и ode23 изо всех сил пытаются интегрировать уравнение. Кроме того, используйте odeset, чтобы передать в постоянном якобиане
и включить отображение статистики решателя.
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.941491 seconds.
title('ode15s')
subplot(2,2,2)
tic, ode23s(@(t,y) -lambda*y, tspan, y0, opts), toc63 successful steps 0 failed attempts 191 function evaluations 0 partial derivatives 63 LU decompositions 189 solutions of linear systems Elapsed time is 0.483929 seconds.
title('ode23s')
subplot(2,2,3)
tic, ode23t(@(t,y) -lambda*y, tspan, y0, opts), toc95 successful steps 0 failed attempts 125 function evaluations 0 partial derivatives 28 LU decompositions 123 solutions of linear systems Elapsed time is 0.632605 seconds.
title('ode23t')
subplot(2,2,4)
tic, ode23tb(@(t,y) -lambda*y, tspan, y0, opts), toc71 successful steps 0 failed attempts 167 function evaluations 0 partial derivatives 23 LU decompositions 236 solutions of linear systems Elapsed time is 0.357163 seconds.
title('ode23tb')
Жесткие решатели, которые все выполняют хорошо, но ode23s завершает интегрирование с наименьшим количеством шагов и запускает самое быстрое для этой конкретной проблемы. Поскольку постоянный якобиан задан, ни один из решателей не должен вычислять частные производные, чтобы вычислить решение. Определение якобиана приносит пользу ode23s больше всего, поскольку это обычно оценивает якобиан на каждом шаге.
Для общих жестких проблем производительность жестких решателей отличается в зависимости от формата проблемы и заданных опций. Обеспечение якобиевской матрицы или шаблона разреженности всегда повышает эффективность решателя для жестких проблем. Но поскольку жесткие решатели используют якобиан по-другому, улучшение может значительно отличаться. В сущности, если система уравнений является очень большой или должна много раз решаться, то стоит исследовать производительность различных решателей, чтобы минимизировать время выполнения.
ode23tb является реализацией TR-BDF2, неявной формулы Рунге-Кутта с шагом метода трапеций как его первая стадия и формула дифференцирования назад порядка два как его второй этап. Конструкцией та же матрица итерации используется в оценке обоих этапов. Как ode23s и ode23t, этот решатель может быть более эффективным, чем ode15s для проблем с грубыми допусками [1], [2].
[1] Банк, Р. E. W. C. Корэн младший W. Fichtner, Э. Гросс, Д. Роуз и Р. Смит, “Переходное Моделирование Кремниевых Устройств и Схем”, Сделка IEEE CAD, 4 (1985), стр 436–451.
[2] Шемпин, L. F. и М. E. Осия, “Анализ и реализация TR-BDF2”, прикладная числовая математика 20, 1996.