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

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

Дифференциальные алгебраические уравнения (DAE) — это дифференциальные уравнения, на эволюцию которых накладываются уравнения связей. В этом руководстве содержатся сведения о функциональных возможностях для решения дифференциальных алгебраических уравнений (DAE). Другие руководства можно найти в SciMLTutorials.jl.

Здесь предполагается, что вы ознакомились с руководством об обыкновенных дифференциальных уравнениях.

Дифференциально-алгебраические уравнения (DAE) матрицы масс

Вместо простого определения ODE в виде можно выразить дифференциальное уравнение в форме с матрицей масс:

где является матрицей масс. Решим уравнение Робертсона. В предыдущих руководствах мы записывали это уравнение в следующем виде:

Но его можно написать с помощью отношения сохранения:

В таком виде его можно записать как ODE матрицы масс, где  — сингуляр (это еще одна форма дифференциально-алгебраического уравнения (DAE)). Здесь последняя строка M просто нулевая. Мы можем реализовать эту форму в виде:

using DifferentialEquations
using Plots
function rober(du, u, p, t)
    y₁, y₂, y₃ = u
    k₁, k₂, k₃ = p
    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
    du[2] = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2
    du[3] = y₁ + y₂ + y₃ - 1
    nothing
end
M = [1.0 0 0
     0 1.0 0
     0 0 0]
f = ODEFunction(rober, mass_matrix = M)
prob_mm = ODEProblem(f, [1.0, 0.0, 0.0], (0.0, 1e5), (0.04, 3e7, 1e4))
sol = solve(prob_mm, Rodas5(), reltol = 1e-8, abstol = 1e-8)

plot(sol, xscale = :log10, tspan = (1e-6, 1e5), layout = (3, 1))

Если матрица масс сингулярна, т. е. система является DAE, вам необходимо убедиться, что выбран решатель, который совместим с DAE.

Неявно определенные дифференциально-алгебраические уравнения

В этом примере мы будем решать уравнение Робертсона в его неявной форме:

Это уравнение представляет собой DAE вида:

которое также известно как дифференциальное уравнение связи, где g — уравнение связи. Модель Робертсона можно записать в следующем виде:

₁

с начальными условиями , , , , и .

Процесс работы с DAE такой же, как и с другими типами уравнений, где все, что нужно знать, — это как определить задачу. Задача DAEProblem указывается путем определения обновления f(out,du,u,p,t) на месте, которое в качестве вывода использует значения для изменения out. Чтобы преобразовать ее в DAE, перенесем все переменные на одну сторону. Итак, функция определяется следующим образом:

function f2(out, du, u, p, t)
    out[1] = -0.04u[1] + 1e4 * u[2] * u[3] - du[1]
    out[2] = +0.04u[1] - 3e7 * u[2]^2 - 1e4 * u[2] * u[3] - du[2]
    out[3] = u[1] + u[2] + u[3] - 1.0
end
f2 (generic function with 1 method)

с начальными условиями

u₀ = [1.0, 0, 0]
du₀ = [-0.04, 0.04, 0.0]
tspan = (0.0, 100000.0)
(0.0, 100000.0)

и создаем задачу DAEProblem:

differential_vars = [true, true, false]
prob = DAEProblem(f2, du₀, u₀, tspan, differential_vars = differential_vars)
DAEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100000.0)
u0: 3-element Vector{Float64}:
 1.0
 0.0
 0.0
du0: 3-element Vector{Float64}:
 -0.04
  0.04
  0.0

differential_vars — это параметр, указывающий, какие из переменных являются дифференциальными, т. е. не чисто алгебраическими (что означает, что их производная появляется в уравнениях невязки). Это необходимо для того, чтобы алгоритм мог находить согласованные начальные условия. Обратите внимание, что первые две переменные определяются их изменениями, а последняя просто определяется уравнением сохранения. Поэтому мы используем differential_vars = [true,true,false].

Как и в других задачах DifferentialEquations, далее следуют команды для решения и построения графика. Здесь мы будем использовать решатель IDA из Sundials:

using Sundials
sol = solve(prob, IDA())
retcode: Success
Interpolation: 3rd order Hermite
t: 145-element Vector{Float64}:
      0.0
      2.1650624290919708e-5
      4.3301248581839416e-5
      8.660249716367883e-5
      0.00017320499432735766
      0.00034640998865471533
      0.000519614982982073
      0.0008660249716367883
      0.0012124349602915037
      0.0015588784138410544
      ⋮
  65892.74549634113
  70410.93208701594
  74929.11867769076
  79447.30526836557
  83965.49185904038
  88483.6784497152
  93001.86504039001
  97520.05163106482
 100000.0
u: 145-element Vector{Vector{Float64}}:
 [1.0, 0.0, 0.0]
 [0.9999991339757784, 8.655379472200322e-7, 4.862743269778565e-10]
 [0.9999982679523074, 1.7296176426658564e-6, 2.4300498836067054e-9]
 [0.9999965359071958, 3.4509400704738393e-6, 1.3152733713218404e-8]
 [0.9999930718263931, 6.840298689437493e-6, 8.787491743885528e-8]
 [0.9999861436803209, 1.3267575081842676e-5, 5.887445972702439e-7]
 [0.9999792156222534, 1.8881934932916107e-5, 1.9024428136186574e-6]
 [0.999965359954415, 2.7374341605089064e-5, 7.265703979760303e-6]
 [0.9999515054512594, 3.223985657676216e-5, 1.6254692163801515e-5]
 [0.9999376511365062, 3.4736651221263134e-5, 2.7612212272492775e-5]
 ⋮
 [0.025758317807329522, 1.0572300963755252e-7, 0.9742415764696608]
 [0.02433027269880141, 9.971731290439181e-8, 0.9756696275838858]
 [0.023046874388090485, 9.433453306501558e-8, 0.9769530312773764]
 [0.021898964410629017, 8.95329107641391e-8, 0.9781009460564604]
 [0.020870482976664065, 8.523980107043642e-8, 0.979129431783535]
 [0.01993615164271255, 8.134667995102116e-8, 0.9800637670106075]
 [0.019084057376200194, 7.780380148123644e-8, 0.9809158648199984]
 [0.018300030088126588, 7.45472133324717e-8, 0.9816998953646601]
 [0.017894798357413172, 7.286692273037954e-8, 0.9821051287756642]

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

plot(sol, xscale = :log10, tspan = (1e-6, 1e5), layout = (3, 1))