Решение дифференциальных алгебраических уравнений (ДАУ)

Что такое дифференциальное алгебраическое уравнение?

Дифференциальные алгебраические уравнения являются типом дифференциального уравнения, где одна или несколько производных зависимых переменных не присутствуют в уравнениях. Переменные, которые появляются в уравнениях без их производной, называются алгебраическими, а наличие алгебраических переменных означает, что вы не можете записать уравнения в явной форме $y' = f (t, y)$ . Вместо этого можно решить ДАУ с этими формами:

ode15s и ode23t решатели могут решить индекс -1 линейно неявные задачи с сингулярной большой матрицы $M (t, y) y' = f (t, y)$ , включая полу-явные ДАУ вида
$\begin{matrix} y' = f (t, y, z) \\ 0 = g (t, y, z) . \end{matrix}$
В этом виде наличие алгебраических переменных приводит к большой матрице, поскольку на основной диагонали есть одно или несколько нулей.
$M y' = (\begin{matrix} y'_{1} & 0 & \dots & 0 \\ 0 & y'_{2} & 0 & ⋮ \\ ⋮ & 0 & ⋱ & 0 \\ 0 & \dots & 0 & 0 \end{matrix}) .$
По умолчанию решатели автоматически проверяют особенность большой матрицы, чтобы обнаружить системы ДАУ. Если вы знаете о особенности раньше времени, то можете задать MassSingular опция odeset на 'yes'. С помощью ДАУ можно также предоставить решателю предположение о начальных условиях для $y'_{0}$ использование InitialSlope свойство odeset. Это в дополнение к определению обычных начальных условий для $y_{0}$ в вызове решателя.
ode15i решатель может решить более общие ДАУ в полностью неявной форме
$f (t, y, y') = 0 .$
В полностью неявной форме наличие алгебраических переменных приводит к сингулярной якобиевой матрице. Это потому, что, по крайней мере, один из столбцов в матрице гарантированно содержит все нули, поскольку производная этой переменной не отображается в уравнениях.
$J = \partial f / \partial y' = (\begin{matrix} \frac{\partial f_{1}}{\partial y'_{1}} & \dots & \frac{\partial f_{1}}{\partial y'_{n}} \\ ⋮ & ⋱ & ⋮ \\ \frac{\partial f_{m}}{\partial y'_{1}} & \dots & \frac{\partial f_{m}}{\partial y'_{n}} \end{matrix})$
ode15i решатель требует, чтобы вы задали начальные условия для обоих $y'_{0}$ и $y_{0}$ . Кроме того, в отличие от других решателей ОДУ, ode15i требуется функция, кодирующая уравнения, чтобы принять дополнительный вход: odefun(t,y,yp).

ДАУ возникают в самых разных системах, потому что законы о физической сохранности часто имеют такие формы, как $x + y + z = 0$ . Если x, x', y, и y' определяются явно в уравнениях, тогда этого уравнения сохранения достаточно, чтобы решить для z без выражения для z'.

Согласованные начальные условия

Когда вы решаете ДАУ, можно задать начальные условия для обоих $y'_{0}$ и $y_{0}$ . ode15i решатель требует, чтобы оба начальных условия были заданы как входные параметры. Для ode15s и ode23t решатели, начальное условие для $y'_{0}$ необязательно (но может быть задано с помощью InitialSlope опция odeset). В обоих случаях возможно, что заданные вами начальные условия не согласуются с уравнениями, которые вы пытаетесь решить. Начальные условия, которые конфликтуют друг с другом, называются несогласованными. Обработка начальных условий варьируется в зависимости от решателя:

ode15s и ode23t - Если вы не задаете начальное условие для $y'_{0}$ , затем решатель автоматически вычисляет согласованные начальные условия на основе начального условия, которое вы обеспечиваете $y_{0}$ . Если вы задаете несогласованное начальное условие для $y'_{0}$ затем решатель обрабатывает значения как догадки, пытается вычислить допустимые значения, близкие к догадкам, и продолжает решать задачу.
ode15i - Начальные условия, которые вы подаете на решатель, должны быть последовательными, и ode15i не проверяет представленные значения на непротиворечивость. Функция помощника decic вычисляет допустимые начальные условия для этой цели.

Дифференциальный индекс

ДАУ характеризуются дифференциальным индексом, который является мерой их особенности. Путем дифференцирования уравнений можно исключить алгебраические переменные, и если вы делаете это достаточно раз, то уравнения принимают форму системы явных ОДУ. Дифференциальный индекс системы ДАУ - это количество производных, которые вы должны взять, чтобы выразить систему как эквивалентную системе явных ОДУ. Таким образом, ОДУ имеют дифференциальный индекс 0.

Примером индекс-1 ДАУ является

$y (t) = k (t) .$

Для этого уравнения можно взять одну производную, чтобы получить явную форму ОДУ

$y' = k' (t) .$

Примером индекс-2 ДАУ является

$\begin{array}{l} y'_{1} = y_{2} \\ 0 = k (t) - y_{1} . \end{array}$

Эти уравнения требуют, чтобы две производные были переписаны в явной форме ОДУ

$\begin{array}{l} y'_{1} = k' (t) \\ y'_{2} = k'' (t) . \end{array}$

ode15s и ode23t решатели решают только ДАУ с индексом 1. Если индекс ваших уравнений 2 или выше, то нужно переписать уравнения как эквивалентную систему ДАУ с индексом-1. Всегда можно взять производные и переписать систему ДАУ как эквивалентную системе индекс-1 ДАУ. Имейте в виду, что если вы заменяете алгебраические уравнения на их производные, то, возможно, вы удалили некоторые ограничения. Если уравнения больше не включают исходные ограничения, то численное решение может дрейфовать.

Если у вас есть Toolbox™ Symbolic Math, то смотрите Решение Дифференциальных Алгебраических Уравнений (ДАУ) (Symbolic Math Toolbox) для получения дополнительной информации.

Навязывание неотрицательности

Большинство опций в odeset работать, как и ожидалось, с решателями ДАУ ode15s, ode23t, и ode15i. Однако одним из заметных исключений является использование NonNegative опция. The NonNegative опция не поддерживает неявные решатели (ode15s, ode23t, ode23tb) применяется к задачам с большой матрицей. Поэтому вы не можете использовать эту опцию, чтобы наложить ограничения неотрицательности на задачу ДАУ, которая обязательно имеет сингулярную большую матрицу. Для получения дополнительной информации см. раздел [1].

Решите задачу Робертсона как полу-явные дифференциальные алгебраические уравнения (ДАУ)

Открыть скрипт

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

$\begin{array}{cl} y'_1 &= -0.04y_1 + 10^4 y_2y_3\\ y'_2 &= 0.04y_1 -
10^4 y_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} y'_1 &= -0.04y_1 + 10^4 y_2y_3\\ 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}$

Дифференциальный индекс этой системы равен 1, поскольку $y_3$ для превращения этой системы в систему ОДУ требуется только одна производная. Поэтому никаких дальнейших преобразований перед решением системы не требуется.

Функция robertsdae кодирует эту систему ДАУ. Сохраните robertsdae.m в текущей папке, чтобы запустить пример.

function out = robertsdae(t,y)
out = [-0.04*y(1) + 1e4*y(2).*y(3)
   0.04*y(1) - 1e4*y(2).*y(3) - 3e7*y(2).^2
   y(1) + y(2) + y(3) - 1 ];

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

Решить систему ДАУ можно используя ode15s. Допустимые начальные условия для y0 очевидны на основе закона о сохранении. Использование odeset чтобы задать опции:

Используйте постоянную большую матрицу, чтобы представлять левую сторону системы уравнений.

$\left( \begin{array}{c} y'_1\\ y'_2\\ 0 \end{array} \right) = M y'
\rightarrow M = \left( \begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 &
0 \end{array} \right)$

Установите допуск относительной погрешности равным 1e-4.
Используйте абсолютную погрешность 1e-10 для второго компонента решения, поскольку шкала существенно изменяется по сравнению с другими компонентами.
Выйдите из 'MassSingular' опция на своём значении по умолчанию 'maybe' для проверки автоматического обнаружения ДАУ.

y0 = [1; 0; 0];
tspan = [0 4*logspace(-6,6)];
M = [1 0 0; 0 1 0; 0 0 0];
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6]);
[t,y] = ode15s(@robertsdae,tspan,y0,options);

Постройте график решения.

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

Ссылки

[1] шемпин, L.F., S. Thompson, J.A. Kierzenka, and G.D. Byrne. «Неотрицательные решения ОДУ». Прикладная математика и расчеты. Том 170, 2005, стр. 556-569.

См. также

ode15i | ode15s | ode23t | odeset

Внешние веб-сайты

Решение Индекса -1 ДАУ в MATLAB и Simulink

Документация