exponenta event banner

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')

Figure contains an axes. The axes contains 2 objects of type line. These objects represent ode15i, exact.

Этот пример переформулирует систему ОДУ как полностью неявную систему дифференциальных алгебраических уравнений (DAE). Проблема Робертсона, закодированная 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.$$

Проблему можно переписать как систему DAE, используя закон сохранения для определения состояния. $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 '1 y2 = 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 дифференциальными алгебраическими уравнениями (DAE). Вспомогательная функция decic вычисляет согласованные начальные условия, которые подходят для использования с ode15i [1].

Ссылки

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

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