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;

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

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

f(t0,y,y)=0.

Поскольку возможно поставить несогласованные начальные условия, и ode15i не проверяет согласованность, рекомендуется использовать функцию helper 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')

Figure contains an axes. The axes contains 2 objects of type line. These objects represent 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 должен принять три входов для 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]. The 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] Lawrence F. Shampine, «Solving 0 = F (t, y (t), y ′ (t)) в MATLAB», Journal of Numerical Mathematics, Vol.10, No.4, 2002, pp. 291-310.

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