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

Что такое Дифференциальное Алгебраическое Уравнение?

Дифференциальные алгебраические уравнения являются типом дифференциального уравнения, где одна или несколько производных зависимых переменных не присутствуют в уравнениях. Переменные, которые появляются в уравнениях без их производной, называются алгебраические, и присутствие алгебраических переменных означает, что вы не можете записать уравнения в явной форме $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$ . Если xXY, и 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. Следует иметь в виду, что, если вы заменяете алгебраические уравнения на их производные, затем вы можете удалить некоторые ограничения. Если уравнения больше не включают исходные ограничения, то числовое решение может дрейфовать.

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

Наложение неотрицательности

Большинство опций в odeset работайте как ожидалось с решателями ДАУ ode15sode23t, и ode15i. Однако одно существенное исключение с использованием NonNegative опция. NonNegative опция не поддерживает неявные решатели (ode15sode23tode23tb) примененный проблемы с большой матрицей. Поэтому вы не можете использовать эту опцию, чтобы наложить ограничения неотрицательности на проблему ДАУ, которая обязательно имеет сингулярную большую матрицу. Для получения дополнительной информации см. [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., С. Томпсон, Дж.А. Кирженка и Г.Д. Бирн. "Неотрицательные решения ОДУ". Прикладная математика и Расчет. Издание 170, 2005, стр 556-569.

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

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

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