24.2 Дифференциально-алгебраические уравнения ¶

Функция daspk может быть использована для решения ДАУ вида

0 = f (x-dot, x, t),    x(t=0) = x_0, x-dot(t=0) = x-dot_0

где x-dot — производная от x. Уравнение решается с помощью решателя ДАУ Petzold’s DASPK.

: [x, xdot, istate, msg] = daspk (fcn, x_0, xdot_0, t, t_crit) ¶

Решить систему дифференциально-алгебраических уравнений.

daspk решает систему уравнений

0 = f (x, xdot, t)

с

x(t_0) = x_0, xdot(t_0) = xdot_0

Решение возвращается в матрицах x и xdot, где каждая строка в результатах соответствует одному из элементов в векторе t. Первый элемент t должен быть t_0 и соответствовать начальному состоянию системы x_0 и его производной xdot_0, так что первая строка выходного x будет x_0, а первая строка выходного xdot — xdot_0.

Первый аргумент, fcn, — строка, встроенная функция или дескриптор функции, которые задают функцию f для вычисления вектора остатков для заданной системы уравнений. Она должна иметь вид

res = f (x, xdot, t)

в котором x, xdot и res — векторы, а t — скаляр.

Если fcn — массив строк из двух элементов или двумерный массив ячеек строк, встроенных функций или дескрипторов функций, то первый элемент задаёт функцию f, описанную выше, а второй — функцию для вычисления модифицированной матрицы Якоби

df       df
jac = -- + c ------
      dx     d xdot

Функция модифицированной матрицы Якоби должна иметь вид

jac = j (x, xdot, t, c)

Второй и третий аргументы к daspk задают начальное состояние состояний и их производных, а четвёртый аргумент — вектор выходных моментов времени, в которых требуется решение, включая момент, соответствующий начальному условию.

Набор начальных состояний и производных не обязательно должен быть согласован. Если они не согласованы, вы должны использовать функцию daspk_options для предоставления дополнительной информации, чтобы daspk мог вычислить согласованную начальную точку.

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

После успешного вычисления значение istate будет больше нуля (согласуется с фортрановской версией DASPK).

Если вычисление не прошло успешно, значение istate будет меньше нуля, и msg будет содержать дополнительную информацию.

Вы можете использовать функцию daspk_options для установки необязательных параметров для daspk.

См. также: dassl.

: daspk_options () ¶

: val = daspk_options (opt) ¶

: daspk_options (opt, val) ¶

Запрос или установка параметров для функции daspk.

При вызове без аргументов отображаются имена всех доступных параметров и их текущие значения.

При передаче одного аргумента возвращается значение параметра opt.

При передаче двух аргументов daspk_options устанавливает параметр opt в значение val.

Параметры включают

"absolute tolerance"

Абсолюльная погрешность. Может быть вектором или скаляром. Если вектором, то он должен соответствовать размеру вектора состояния, а относительная погрешность также должна быть вектором той же длины.

"relative tolerance"

Относительная погрешность. Может быть вектором или скаляром. Если вектором, то он должен соответствовать размеру вектора состояния, а абсолютная погрешность также должна быть вектором той же длины.

Локальный тест погрешности, применяемый на каждом шаге интегрирования,

abs (local error in x(i))
       <= rtol(i) * abs (Y(i)) + atol(i)

"compute consistent initial condition"

Обозначая дифференциальные переменные в векторе состояния как ‘Y_d’, а алгебраические переменные как ‘Y_a’, ddaspk может решить одну из двух задач инициализации:

1. При заданном Y_d вычислить Y_a и Y’_d
2. При заданном Y’ вычислить Y.

В любом случае, входными данными являются начальные значения заданных компонентов, а также должны быть предоставлены начальные предположения для неизвестных компонентов. Установите этот параметр в 1, чтобы решить первую задачу, или в 2, чтобы решить вторую (по умолчанию 0, поэтому вы должны предоставить набор начальных условий, которые согласованы).

Если этот параметр установлен на ненулевое значение, вы также должны установить параметр "algebraic variables" для объявления алгебраических переменных в проблеме.

"use initial condition heuristics"

Установите это значение в ненулевое, чтобы использовать варианты эвристики начального условия, описанные ниже.

"initial condition heuristics"

Вектор параметров, которые можно использовать для управления вычислением начального условия.

MXNIT

Максимальное количество итераций Ньютона (по умолчанию 5).

MXNJ

Максимальное количество вычислений Якоби (по умолчанию 6).

MXNH

Максимальное количество значений искусственного параметра шага, которые нужно попробовать, если параметр "compute consistent initial condition" установлен в 1 (по умолчанию 5).

Обратите внимание, что максимальное общее число итераций Ньютона, разрешённое, равно MXNIT*MXNJ*MXNH если параметр "compute consistent initial condition" установлен в 1 и MXNIT*MXNJ если он установлен в 2.

LSOFF

Установите в ненулевое значение, чтобы отключить алгоритм поиска по линии (по умолчанию 0).

STPTOL

Минимальный масштабированный шаг в алгоритме поиска по линии (по умолчанию eps^(2/3)).

EPINIT

Фактор качания в тесте сходимости итерации Ньютона. Тест применяется к вектору остатков, умноженному слева на приближённую матрицу Якоби. Для сходимости взвешенная среднеквадратическая норма этого вектора (масштабированная весовыми коэффициентами) должна быть меньше EPINIT*EPCON, где EPCON = 0,33 — аналогичная константа теста, используемая в шагах по времени. По умолчанию EPINIT = 0,01.

"print initial condition info"

Установите этот параметр в ненулевое значение, чтобы отобразить подробную информацию о расчёте начального состояния (по умолчанию 0).

"exclude algebraic variables from error test"

Установите в ненулевое значение, чтобы исключить алгебраические переменные из теста погрешности. Вы также должны установить параметр "algebraic variables" для объявления алгебраических переменных в проблеме (по умолчанию 0).

"algebraic variables"

Вектор той же длины, что и вектор состояния. ненулевые элементы указывают, что соответствующий элемент вектора состояния является алгебраической переменной (т. е. его производная не появляется явно в наборе уравнений).

Этот параметр необходим для параметров "compute consistent initial condition" и "exclude algebraic variables from error test".

"enforce inequality constraints"

Установите это значение в одно из следующих, чтобы принудительно выполнить неравенства, заданные параметром "inequality constraint types" (по умолчанию 0).

1. Выполнять проверку ограничений только при вычислении начального состояния.
2. Выполнять проверку ограничений во время интегрирования.
3. Выполнять обе опции 1 и 2.

"inequality constraint types"

Вектор той же длины, что и состояние, определяющий тип неравенства. Каждый элемент вектора соответствует элементу состояния и должен быть назначен одному из следующих кодов

-2

Меньше нуля.

-1

Меньше или равно нулю.

0

Без ограничений.

1

Больше или равно нулю.

2

Больше нуля.

Этот параметр действует только если параметр "enforce inequality constraints" ненулевой.

"initial step size"

Дифференциально-алгебраические проблемы иногда могут испытывать серьёзные трудности с масштабированием на первом шаге. Если вам много известно о масштабировании вашей проблемы, вы можете помочь смягчить эту проблему, задав начальный шаг (по умолчанию вычисляется автоматически).

"maximum order"

Ограничьте максимальный порядок метода решения. Этот параметр должен быть между 1 и 5 включительно (по умолчанию 5).

"maximum step size"

Установка максимального шага позволит избежать прохождения очень больших областей (значение по умолчанию не задано).

Octave также включает в себя DASSL, более раннюю версию DASPK, и DASRT, которые могут быть использованы для решения ДАУ с ограничениями (условиями остановки).

: [x, xdot, istate, msg] = dassl (fcn, x_0, xdot_0, t, t_crit) ¶

Решить систему дифференциально-алгебраических уравнений.

dassl решает систему уравнений

0 = f (x, xdot, t)

с

x(t_0) = x_0, xdot(t_0) = xdot_0

Решение возвращается в матрицах x и xdot, где каждая строка в результирующих матрицах соответствует одному из элементов вектора t. Первый элемент t должен быть t_0 и соответствовать начальному состоянию системы x_0 и её производной xdot_0, так что первая строка выходной x — x_0, а первая строка выходной xdot — xdot_0.

Первый аргумент, fcn, — строка, встроенная функция или дескриптор функции, которые называют функцию f для вычисления вектора остатков для системы уравнений. Она должна иметь вид

res = f (x, xdot, t)

где x, xdot и res — векторы, а t — скаляр.

Если fcn — двухэлементный массив строк или двухэлементный массив ячеек строк, встроенных функций или дескрипторов функций, первый элемент называет функцию f, описанную выше, а второй элемент называет функцию для вычисления модифицированной якобиана

df       df
jac = -- + c ------
      dx     d xdot

Функция модифицированного якобиана должна иметь вид

jac = j (x, xdot, t, c)

Второй и третий аргументы для dassl задают начальное состояние состояний и их производных, а четвёртый аргумент задаёт вектор выходных моментов времени, где требуется решение, включая момент времени, соответствующий начальному условию.

Множество начальных состояний и производных не обязательно должны быть согласованными. Однако на практике DASSL не очень хорошо определяет согласованный набор, поэтому лучше убедиться, что начальные значения приводят к тому, что функция равна нулю.

Пятый аргумент необязательный и может использоваться для задания набора моментов времени, до которых решатель DAE не должен интегрировать. Он полезен для предотвращения трудностей с особенностями и точками, где есть разрыв в производной.

После успешного вычисления значение istate будет больше нуля (согласуется с Fortran-версией DASSL).

Если вычисление не прошло успешно, значение istate будет меньше нуля, а msg будет содержать дополнительную информацию.

Вы можете использовать функцию dassl_options для установки необязательных параметров для dassl.

См. также: daspk, dasrt, lsode.

: dassl_options () ¶

: val = dassl_options (opt) ¶

: dassl_options (opt, val) ¶

Запрос или установка параметров для функции dassl.

При вызове без аргументов отображаются имена всех доступных параметров и их текущие значения.

При передаче одного аргумента возвращается значение параметра opt.

При вызове с двумя аргументами, dassl_options устанавливает параметр opt в значение val.

Параметры включают

"absolute tolerance"

Абсолюльная погрешность. Может быть вектором или скаляром. Если вектором, то он должен соответствовать размерности вектора состояния, а относительная погрешность также должна быть вектором той же длины.

"relative tolerance"

Относительная погрешность. Может быть вектором или скаляром. Если вектором, то он должен соответствовать размерности вектора состояния, а абсолютная погрешность также должна быть вектором той же длины.

Локальный тест ошибки, применяемый на каждом шаге интегрирования,

abs (local error in x(i))
       <= rtol(i) * abs (Y(i)) + atol(i)

"compute consistent initial condition"

Если не равно нулю, dassl попытается вычислить согласованный набор начальных условий. Как правило, это ненадёжно, поэтому лучше задавать согласованный набор и оставлять этот параметр равным нулю.

"enforce nonnegativity constraints"

Если известно, что решения ваших уравнений всегда будут неотрицательными, это может помочь установить этот параметр в ненулевое значение. Однако, вероятно, лучше сначала попробовать оставить этот параметр равным нулю и установить его в ненулевое значение только в том случае, если это не работает очень хорошо.

"initial step size"

Дифференциально-алгебраические задачи иногда могут испытывать серьёзные проблемы со масштабированием на первом шаге. Если вы хорошо знаете масштабирование вашей задачи, вы можете помочь решить эту проблему, задав начальный шаг.

"maximum order"

Ограничить максимальный порядок метода решения. Этот параметр должен быть в диапазоне от 1 до 5 включительно.

"maximum step size"

Установление максимального шага позволит избежать пропуска очень больших областей (по умолчанию не указано).

"step limit"

Максимальное количество шагов интегрирования, которые нужно предпринять при одном вызове основной Fortran-функции.

: [x, xdot, t_out, istat, msg] = dasrt (fcn, g, x_0, xdot_0, t) ¶

: … = dasrt (fcn, g, x_0, xdot_0, t, t_crit) ¶

: … = dasrt (fcn, x_0, xdot_0, t) ¶

: … = dasrt (fcn, x_0, xdot_0, t, t_crit) ¶

Решает систему дифференциально-алгебраических уравнений.

dasrt решает систему уравнений

0 = f (x, xdot, t)

с

x(t_0) = x_0, xdot(t_0) = xdot_0

с функциональными критериями остановки (решение корней).

Решение возвращается в матрицах x и xdot, где каждая строка в результирующих матрицах соответствует одному из элементов в векторе t_out. Первый элемент t должен быть t_0 и соответствовать начальному состоянию системы x_0 и её производной xdot_0, так что первая строка результата x — x_0, а первая строка результата xdot — xdot_0.

Вектор t задаёт верхний предел длины интегрирования. Если условие остановки выполнено, вектор t_out будет короче t, и последний элемент t_out будет точкой, в которой было выполнено условие остановки, и может не соответствовать какому-либо элементу вектора t.

Первый аргумент, fcn, — строка, встроенная функция или дескриптор функции, которая называет функцию f для вычисления вектора остатков для набора уравнений. Она должна иметь вид

res = f (x, xdot, t)

в котором x, xdot и res — векторы, а t — скаляр.

Если fcn — строковый массив или ячейный массив из двух строк, встроенных функций или дескрипторов функций, первый элемент называет функцию f, описанную выше, а второй элемент называет функцию для вычисления модифицированной якобиана

df       df
jac = -- + c ------
      dx     d xdot

Функция модифицированного якобиана должна иметь вид

jac = j (x, xdot, t, c)

Необязательный второй аргумент называет функцию, которая определяет функции ограничений, корни которых желательно найти во время интегрирования. Эта функция должна иметь вид

g_out = g (x, t)

и возвращать вектор значений функций ограничений. Если значение любой из функций ограничений меняет знак, DASRT попытается остановить интегрирование в точке изменения знака.

Если имя функции ограничения опущено, dasrt решает ту же проблему, что и daspk или dassl.

Обратите внимание, что из-за погрешностей в функциях ограничений, связанных с округлением и ошибкой интегрирования, DASRT может возвращать ложные корни или возвращать один и тот же корень в двух или более близких значениях T. Если предполагается наличие таких ложных корней, пользователь должен рассмотреть более низкие значения толерантностей или более высокую точность при вычислении функций ограничений.

Если корень какой-либо функции ограничения определяет конец задачи, входные данные для DASRT должны тем не менее позволять интегрирование в точку немного за этим корнем, чтобы DASRT мог найти корень интерполяцией.

Третий и четвёртый аргументы dasrt задают начальное состояние состояний и их производных, а четвёртый аргумент задаёт вектор времен выхода, в которых требуется решение, включая время, соответствующее начальному состоянию.

Набор начальных состояний и производных не обязательно должен быть согласованным. Однако на практике DASSL не очень хорошо определяет согласованный набор для вас, поэтому лучше убедиться, что начальные значения приводят к нулевому значению функции.

Шестой аргумент необязательный и может быть использован для задания набора времен, до которых решатель DAE не должен интегрировать. Это полезно для избежания проблем с особенностями и точками, где есть разрыв в производной.

После успешного вычисления значение istate будет больше нуля (согласуясь с Fortran-версией DASSL).

Если вычисление не было успешным, значение istate будет меньше нуля, а msg будет содержать дополнительную информацию.

Вы можете использовать функцию dasrt_options для настройки необязательных параметров для dasrt.

См. также: dasrt_options, daspk, dasrt, lsode.

: dasrt_options () ¶

: val = dasrt_options (opt) ¶

: dasrt_005foptions (opt, val) ¶

Запрос или установка параметров для функции dasrt.

При вызове без аргументов отображаются имена всех доступных параметров и их текущие значения.

При вызове с одним аргументом возвращается значение параметра opt.

При вызове с двумя аргументами dasrt_options устанавливает параметр opt в значение val.

Параметры включают

"absolute tolerance"

Абсолютная точность. Может быть как вектором, так и скаляром. Если вектором, то он должен соответствовать размеру вектора состояния, и относительная точность также должна быть вектором той же длины.

"relative tolerance"

Относительная точность. Может быть как вектором, так и скаляром. Если вектором, то он должен соответствовать размеру вектора состояния, и абсолютная точность также должна быть вектором той же длины.

Локальный тест на ошибку, применяемый на каждом шаге интегрирования, —

abs (local error in x(i)) <= ...
      rtol(i) * abs (Y(i)) + atol(i)

"initial step size"

Дифференциально-алгебраические задачи иногда могут испытывать значительные трудности со масштабированием на первом шаге. Если вы хорошо знаете масштабирование вашей задачи, вы можете помочь решить эту проблему, указав начальный шаг.

"maximum order"

Ограничение максимального порядка метода решения. Этот параметр должен быть от 1 до 5 включительно.

"maximum step size"

Установка максимального шага позволит избежать прохождения очень больших областей.

"step limit"

Максимальное число шагов интегрирования, которые необходимо попытаться выполнить за один вызов кода Fortran.

См. K. E. Brenan и др., Численное решение начальных задач в дифференциально-алгебраических уравнениях, North-Holland (1989), DOI: https://doi.org/10.1137/1.9781611971224, для получения более подробной информации о реализации DASSL.