Engee documentation

Differential Algebraic Equations

Differential Algebraic Equations (DAEs) are differential equations which have constraint equations on their evolution. This tutorial will introduce you to the functionality for solving differential algebraic equations (DAEs). Other introductions can be found by checking out SciMLTutorials.jl.

This tutorial assumes you have read the Ordinary Differential Equations tutorial.

Mass-Matrix Differential-Algebraic Equations (DAEs)

Instead of just defining an ODE as , it can be common to express the differential equation in the form with a mass matrix:

where is known as the mass matrix. Let’s solve the Robertson equation. In previous tutorials, we wrote this equation as:

But we can instead write this with a conservation relation:

In this form, we can write this as a mass matrix ODE where is singular (this is another form of a differential-algebraic equation (DAE)). Here, the last row of M is just zero. We can implement this form as:

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

If your mass matrix is singular, i.e. your system is a DAE, then you need to make sure you choose a solver that is compatible with DAEs

Implicitly-Defined Differential-Algebraic Equations (DAEs)

In this example, we will solve the Robertson equation in its implicit form:

This equation is a DAE of the form:

which is also known as a constrained differential equation, where g is the constraint equation. The Robertson model can be written in the form:

₁

with initial conditions , , , , , and .

The workflow for DAEs is the same as for the other types of equations, where all you need to know is how to define the problem. A DAEProblem is specified by defining an in-place update f(out,du,u,p,t) which uses the values to mutate out as the output. To makes this into a DAE, we move all the variables to one side. Thus, we can define the function:

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)

with initial conditions

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

and make the 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 is an option which states which of the variables are differential, i.e. not purely algebraic (which means that their derivative shows up in the residual equations). This is required for the algorithm to be able to find consistent initial conditions. Notice that the first two variables are determined by their changes, but the last is simply determined by the conservation equation. Thus, we use differential_vars = [true,true,false].

As with the other DifferentialEquations problems, the commands are then to solve and plot. Here we will use the IDA solver from 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]

In order to clearly see all the features of this solution, it should be plotted on a logarithmic scale. We’ll also plot each on a different subplot, to allow scaling the y-axis appropriately.

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