A simple but powerful header-only C++ DAE (Differential Algebraic Equation) system solver
A simple but powerful header-only C++ solver for systems of Differential-Algebraic Equations (DAE).
NOTE: dae-cpp
has been redesigned and there were breaking changes between v1.x
and v2.x
. If your project still relies on the old dae-cpp
(v1.x
), it is archived in the legacy branch. For the new version (v2.x
), see Documentation and the notes below.
dae-cpp
dae-cpp
is a cross-platform, header-only C++17 library for solving stiff systems of DAEs (an initial value problem). DAE systems can contain both differential and algebraic equations and can be written in the following matrix-vector form:
$$\mathbf{M}(t) \frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t} = \mathbf{f}(\mathbf{x}, t),$$
to be solved in the interval $t \in [0, t_\mathrm{end}]
$ with the initial condition $\mathbf{x}\rvert_{t=0} = \mathbf{x}_0
$. Here $\mathbf{M}(t)
$ is the mass matrix (can depend on time), $\mathbf{x}(t)
$ is the state vector, and $\mathbf{f}(\mathbf{x}, t)
$ is the (nonlinear) vector function of the state vector $\mathbf{x}
$ and time $t$.
The DAE solver uses implicit Backward Differentiation Formulae (BDF) of orders I-IV with adaptive time stepping. Every time step, the BDF integrator reduces the original DAE system to a system of nonlinear equations, which is solved using iterative Quasi-Newton root-finding algorithm. The Quasi-Newton method reduces the problem further down to a system of linear equations, which is solved using Eigen, a versatile and fast C++ template library for linear algebra. Eigen's sparse solver performs two steps: factorization (decomposition) of the Jacobian matrix and the linear system solving itself. This gives us the numerical solution of the entire DAE system at the current time step. Finally, depending on the convergence rate of the Quasi-Newton method, variability of the solution, and user-defined accuracy, the DAE solver adjusts the time step size and initiates a new iteration in time.
This library is header only, no need to install, just copy dae-cpp
, Eigen
, and autodiff
folders into your project.
Examples and tests can be compiled using CMake (see Testing).
If you already have cloned the project without --recurse-submodules
option, you can initialize and update googletest
submodule by running
git submodule update --init
Then build and run the tests:
mkdir build && cd build
cmake ..
make
ctest
Consider the following (trivial) DAE system as a quick example:
\left\{
\begin{alignedat}{3}
\dot x & = y, \\
y & = \cos(t),
\end{alignedat}
\right.
with the initial condition:
\left\{
\begin{alignedat}{3}
x\rvert_{t=0} & = 0, \\
y\rvert_{t=0} & = 1.
\end{alignedat}
\right.
This system contains one simple differential equation and one algebraic equation. The analytic solution is the following:
\left\{
\begin{alignedat}{3}
x(t) & = \sin(t), \\
y(t) & = \cos(t).
\end{alignedat}
\right.
Below is a simplified procedure of defining and solving the DAE system using dae-cpp
.
dae-cpp
header into the project#include <dae-cpp/solver.hpp>
Optionally, add daecpp
namespace:
using namespace daecpp;
The mass matrix contains only one non-zero element:
$$ \mathbf{M} = \begin{vmatrix} 1 & 0 \ 0 & 0 \end{vmatrix}. $$
struct MyMassMatrix
{
void operator()(sparse_matrix &M, const double t)
{
M(0, 0, 1.0); // Row 0, column 0, non-zero element 1.0
}
};
struct MyRHS
{
void operator()(state_type &f, const state_type &x, const double t)
{
f[0] = x[1]; // y
f[1] = cos(t) - x[1]; // cos(t) - y
}
};
MyMassMatrix mass; // Mass matrix object
MyRHS rhs; // Vector-function object
System my_system(mass, rhs); // Defines the DAE system object
state_vector x0{0, 1}; // The initial state vector (initial condition)
double t{1.0}; // The integration interval: t = [0, 1.0]
my_system.solve(x0, t); // Solves the system with initial condition `x0` and time `t`
or simply
my_system.solve({0, 1}, 1.0);
Vector of solution vectors x
and vector of the corresponding times t
will be stored in my_system.sol.x
and my_system.sol.t
, respectively.
The entire source code is provided in the Quick Start example.
For more information, refer to the Documentation.
Differentiating the RHS w.r.t. $x$ and $y$ gives the following Jacobian matrix:
$$ \mathbf{J} = \begin{vmatrix} 0 & 1 \ 0 & -1 \end{vmatrix}. $$
This matrix can be defined in the code as
struct MyJacobian
{
void operator()(sparse_matrix &J, const state_vector &x, const double t)
{
J.reserve(2); // Pre-allocates memory for 2 non-zero elements (optional)
J(0, 1, 1.0); // Row 0, column 1, non-zero element 1.0
J(1, 1, -1.0); // Row 1, column 1, non-zero element -1.0
}
};
Then add the user-defined Jacobian to the solve()
method:
my_system.solve(x0, t, MyJacobian());
Defining analytic Jacobian matrix can significantly speed up the computation (especially for big systems).
For example, restrict the maximum time step:
my_system.opt.dt_max = 0.1; // Update `dt_max`
my_system.solve(x0, t, MyJacobian()); // Restart the computation
Thank you for considering contributing to the project! Whether you're an experienced developer or just starting out, your ideas and improvements can make this project even better. No contribution is too small!
Feel free to create a GitHub issue for any questions, suggestions, or feedback you may have.