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

• Header only, no pre-compilation required.
• Uses automatic (algorithmic, exact) differentiation (autodiff package) to compute the Jacobian matrix, if it is not provided by the user.
• Fourth-order variable-step implicit BDF time integrator that preserves accuracy even when the time step rapidly changes.
• A very flexible and customizable variable time stepping algorithm based on the solution stability and variability.
• Mass matrix can be non-static (can depend on time) and it can be singular.
• The library is extremely easy to use. A simple DAE can be set up using just a few lines of code (see Quick Start example below).

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

• For more information about the solver, please refer to the Documentation pages.
• Ready to use examples are given in the examples directory of this repository.
• All notable user-facing changes to this project are documented in the CHANGELOG.

Consider the following (trivial) DAE system as a quick example:

with the initial condition:

This system contains one simple differential equation and one algebraic equation. The analytic solution is the following:

Below is a simplified procedure of defining and solving the DAE system using dae-cpp.

#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!

0. Create a GitHub issue if you want to suggest or discuss your changes.
1. Fork the repository and clone it to your local machine.
2. Create a new branch for your contributions.
3. Make your changes and ensure they adhere to our coding standards.
4. Add tests and examples (if relevant), test your changes thoroughly.
5. Submit a pull request with a clear description of your changes and why they are beneficial.

Feel free to create a GitHub issue for any questions, suggestions, or feedback you may have.

• dae-cpp is licensed under the MIT license.
• autodiff is licensed under the MIT license.
• Eigen is licensed under the MPL2 license.