ode15i

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

Синтаксис

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

Описание

пример

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

пример

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

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

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

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

Примеры

свернуть все

Вычислите сопоставимые начальные условия и решите неявный ОДУ с ode15i.

Уравнение Вайссингера

ty2(y)3-y3(y)2+t(t2+1)y-t2y=0.

Поскольку уравнение находится в типичной форме f(t,y,y)=0, можно использовать функцию ode15i, чтобы решить неявное дифференциальное уравнение.

Уравнение кода

Чтобы закодировать уравнение в форме, подходящей для ode15i, необходимо записать функцию с входными параметрами для t, y, и y это возвращает остаточное значение уравнения. Функциональный @weissinger кодирует это уравнение. Просмотрите файл функции.

type weissinger
function res = weissinger(t,y,yp)
%WEISSINGER  Evaluate the residual of the Weissinger implicit ODE
%
%   See also ODE15I.

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

res = t*y^2 * yp^3 - y^3 * yp^2 + t*(t^2 + 1)*yp - t^2 * y;

Вычислите сопоставимые начальные условия

Решатель ode15i требует сопоставимых начальных условий, то есть, начальные условия, предоставленные решателю, должны удовлетворить

f(t0,y,y)=0.

Поскольку возможно предоставить противоречивые начальные условия, и ode15i не проверяет на непротиворечивость, рекомендуется, чтобы вы использовали функцию помощника decic, чтобы вычислить такие условия. decic содержит некоторые заданные зафиксированные переменные и вычисляет сопоставимые начальные значения для незакрепленных переменных.

В этом случае зафиксируйте начальное значение y(t0)=32 и позвольте decic вычислить сопоставимое начальное значение для производной y(t0), запуск с исходного предположения y(t0)=0.

t0 = 1;
y0 = sqrt(3/2);
yp0 = 0;
[y0,yp0] = decic(@weissinger,t0,y0,1,yp0,0)
y0 = 1.2247
yp0 = 0.8165

Решите уравнение

Используйте сопоставимые начальные условия, возвращенные decic с ode15i, чтобы решить ОДУ по временному интервалу [1 10].

[t,y] = ode15i(@weissinger,[1 10],y0,yp0);

Постройте результаты

Точное решение этого ОДУ

y(t)=t2+12.

Постройте числовое решение y, вычисленный ode15i против аналитического решения ytrue.

ytrue = sqrt(t.^2 + 0.5);
plot(t,y,'*',t,ytrue,'-o')
legend('ode15i', 'exact')

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

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

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

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

Функциональный robertsidae кодирует эту систему ДАУ.

function res = robertsidae(t,y,yp)
res = [yp(1) + 0.04*y(1) - 1e4*y(2)*y(3);
   yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2;
   y(1) + y(2) + y(3) - 1];

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

Установите ошибочные допуски и значение.

options = odeset('RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6], ...
   'Jacobian',{[],[1 0 0; 0 1 0; 0 0 0]});

Используйте decic, чтобы вычислить сопоставимые начальные условия из предположений. Зафиксируйте первые два компонента y0, чтобы получить те же сопоставимые начальные условия, как найдено ode15s в hb1dae.m, который формулирует эту проблему как полуявную систему ДАУ.

y0 = [1; 0; 1e-3];
yp0 = [0; 0; 0];
[y0,yp0] = decic(@robertsidae,0,y0,[1 1 0],yp0,[],options);

Решите систему ДАУ с помощью ode15i.

tspan = [0 4*logspace(-6,6)];
[t,y] = ode15i(@robertsidae,tspan,y0,yp0,options);

Постройте компоненты решения. Поскольку второй компонент решения является небольшим относительно других, умножьте его на 1e4 перед графическим выводом.

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

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

свернуть все

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

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

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

function f = odefun(t,y,yp)
f = yp - y;

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

y'1y2=0y'2+1=0 ,

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

function dy = odefun(t,y,yp)
dy = zeros(2,1);
dy(1) = yp(1)-y(2);
dy(2) = yp(2)+1;

Для получения информации о том, как предоставить дополнительные параметры функциональному 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

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

Начальные условия для y0 и yp0 должны быть сопоставимыми, означая это f(t0,y0,y'0)=0. Используйте функцию decic, чтобы вычислить сопоставимые начальные условия близко к предполагаемым значениям.

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

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

Начальные условия для y0 и yp0 должны быть сопоставимыми, означая это f(t0,y0,y'0)=0. Используйте функцию decic, чтобы вычислить сопоставимые начальные условия близко к предполагаемым значениям.

Типы данных: 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. Значения указывают, какое событие решатель обнаружил.

Советы

  • Обеспечение якобиевской матрицы к ode15i очень важно для надежности и эффективности. Также, если система является большой и разреженной, то обеспечение якобиевского шаблона разреженности также помогает решателю. Или в случае, используйте odeset, чтобы передать в матрицах с помощью Jacobian или в опций JPattern.

Алгоритмы

ode15i является переменным шагом, переменный порядок (VSVO) решатель на основе формул дифференцирования назад (BDF) порядков 1 - 5. ode15i разработан, чтобы использоваться с полностью неявными дифференциальными уравнениями и индексным 1 дифференциалом алгебраические уравнения (ДАУ). decic функции помощника вычисляет сопоставимые начальные условия, которые подходят, чтобы использоваться с ode15i [1].

Ссылки

[1] Лоуренс Ф. Шемпин, “Решая 0 = F (t, y (t), y ′ (t)) в MATLAB”, Журнал Числовой Математики, Vol.10, № 4, 2002, стр 291-310.

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