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, и создание соответствующей функции: Значение, isterminal, direction] = myEventFcnTYyp ). Для получения дополнительной информации смотрите Местоположение События ОДУ.

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
%   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 решать ОДУ по временному интервалу [110].

[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, является классическая тестовая задача для программ, которые решают жесткие ОДУ. Система уравнений

$$\begin{array}{cl} y'_1 &= -0.04y_1 + 10^4 y_2y_3\\ y'_2 &= 0.04y_1 -
10^4 4y_2y_3-(3 \times 10^7)y_2^2\\ y'_3 &= (3 \times
10^7)y_2^2.\end{array}$$

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

$$y'_1 + y'_2 + y'_3 = 0.$$

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

$$y_1 + y_2 + y_3 = 1.$$

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

$$\begin{array}{cl} 0 &= y'_1 +0.04y_1 - 10^4 y_2y_3\\ 0 &= y'_2 -0.04y_1
+ 10^4 y_2y_3+(3 \times 10^7)y_2^2\\ 0 &= y_1 + y_2 + y_3 -
1.\end{array}$$

Функциональный 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.

Установите ошибочные погрешности и значение$\partial f / \partial y'$.

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 должен принять три входных параметров для tY, и 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