C++ Solver for Stiff Nonlinear System#
Here’s an implementation of a first-order backward-Euler solver for the same nonlinear stiff system.
This uses a two-dimensional multidimensional array class, Array defined
in the header array.H. This allows us to do things like:
Array A{{{1, 1, 1},
{-1, 2, 0},
{2, 0, 1}}};
to initialize a 2-d array. We can then get the number of rows and
columns via A.nrows() and A.ncols(). We can also index the
array using A(i,j) to get the i-th row and j-th column.
Here’s the class header array.H:
#ifndef ARRAY_H
#define ARRAY_H
#include <vector>
#include <iostream>
#include <cassert>
// a contiguous 2-d array
// here the data is stored in row-major order in a 1-d memory space
// managed as a vector. We overload () to allow us to index this as
// a(irow, icol)
struct Array {
int _rows;
int _cols;
std::vector<double> _data;
Array (int rows, int cols, double val=0.0)
: _rows{rows},
_cols{cols},
_data(rows * cols, val)
{
// we do the asserts here after the initialization of _data
// in the initialization list, but if the size is zero, we'll
// abort here anyway.
assert (rows > 0);
assert (cols > 0);
}
// constructor that allows us to initialize array
// via an initializer list of rows
Array (std::vector<std::vector<double>> v)
: _rows{static_cast<int>(v.size())},
_cols{static_cast<int>(v[0].size())},
_data(_rows * _cols, 0.0)
{
int idx = 0;
for (int i = 0; i < _rows; ++i) {
assert (static_cast<int>(v[i].size()) == _cols);
for (int j = 0; j < _cols; ++j) {
_data[idx] = v[i][j];
++idx;
}
}
}
// note the "const" after the argument list here -- this means
// that this can be called on a const Array
inline int ncols() const { return _cols;}
inline int nrows() const { return _rows;}
inline double& operator()(int row, int col) {
assert (row >= 0 && row < _rows);
assert (col >= 0 && col < _cols);
return _data[row*_cols + col];
}
inline const double& operator()(int row, int col) const {
assert (row >= 0 && row < _rows);
assert (col >= 0 && col < _cols);
return _data[row*_cols + col];
}
};
// the << operator is not part of the of the class, so it is not a member
std::ostream& operator<< (std::ostream& os, const Array& a) {
for (int row = 0; row < a.nrows(); ++row) {
for (int col = 0; col < a.ncols(); ++col) {
os << a(row, col) << " ";
}
os << std::endl;
}
return os;
}
#endif
The Gaussian elimination function is in the header gauss.H:
#ifndef GAUSS_H
#define GAUSS_H
#include <array.H>
#include <vector>
#include <limits>
#include <iostream>
#include <iomanip>
inline
std::vector<double> gauss_elim(Array& A, std::vector<double>& b) {
// perform gaussian elimination with pivoting, solving A x = b.
//
// A is an NxN matrix, x and b are an N-element vectors. Note: A and b are
// changed upon exit to be in upper triangular (row echelon) form
const int N = b.size();
// A is square, with each dimension of length N
assert(A.nrows() == N && A.ncols() == N);
// allocate the solution array
std::vector<double> x(N, 0.0);
// find the scale factors for each row -- this is used when pivoting
std::vector<double> scales(N, 0.0);
for (int irow = 0; irow < N; ++irow) {
double maxval = std::numeric_limits<double>::lowest();
for (int jcol = 0; jcol < N; ++jcol) {
maxval = std::max(maxval, std::abs(A(irow, jcol)));
}
scales[irow] = maxval;
}
// keep track of the number of times we swapped rows
int num_row_swap{};
// main loop over rows
for (int krow = 0; krow < N; ++krow) {
// find the pivot row based on the size of column k -- only consider the
// rows beyond the current row
int row_max{krow};
double col_max = std::numeric_limits<double>::lowest();
for (int kk = krow; kk < N; ++kk) {
double scaled_col_val = std::abs(A(kk, krow)) / scales[kk];
if (scaled_col_val > col_max) {
col_max = scaled_col_val;
row_max = kk;
}
}
// swap the row with the largest scaled element in the current column
// with the current row (pivot) -- do this with b too!
if (row_max != krow) {
// swap row krow with row row_max
for (int jcol = 0; jcol < N; ++jcol) {
double tmp_val = A(krow, jcol);
A(krow, jcol) = A(row_max, jcol);
A(row_max, jcol) = tmp_val;
}
// now swap the same elements in b
double tmp_b = b[krow];
b[krow] = b[row_max];
b[row_max] = tmp_b;
// finally swap the scales
double tmp_scale = scales[krow];
scales[krow] = scales[row_max];
scales[row_max] = tmp_scale;
++num_row_swap;
}
// do the forward-elimination for all rows below the current
for (int irow = krow+1; irow < N; ++irow) {
double coeff = A(irow, krow) / A(krow, krow);
for (int jcol = krow+1; jcol < N; ++jcol) {
A(irow, jcol) += -A(krow, jcol) * coeff;
}
A(irow, krow) = 0.0;
b[irow] += -coeff * b[krow];
}
}
// back-substitution
// last solution is easy
x[N-1] = b[N-1] / A(N-1, N-1);
for (int irow = N-2; irow >= 0; --irow) {
double bsum = b[irow];
for (int jcol = irow+1; jcol < N; ++jcol) {
bsum += -A(irow, jcol) * x[jcol];
}
x[irow] = bsum / A(irow, irow);
}
return x;
}
#endif
and the full solver is in stiff_system.cpp:
#include <cmath>
#include <iostream>
#include <limits>
#include <array.H>
#include <gauss.H>
///
/// given the state Y, return the RHS of our ODE system
///
std::vector<double> rhs(const std::vector<double>& Y) {
std::vector<double> dYdt(3, 0.0);
dYdt[0] = -0.04 * Y[0] + 1.e4 * Y[1] * Y[2];
dYdt[1] = 0.04 * Y[0] - 1.e4 * Y[1] * Y[2] - 3.e7 * Y[1] * Y[1];
dYdt[2] = 3.e7 * Y[1] * Y[1];
return dYdt;
}
///
/// given the state Y, return the Jacobian of the ODE system
///
Array jacobian(const std::vector<double>& Y) {
Array jac(3, 3);
jac(0, 0) = -0.04;
jac(0, 1) = 1.e4 * Y[2];
jac(0, 2) = 1.e4 * Y[1];
jac(1, 0) = 0.04;
jac(1, 1) = -1.e4 * Y[2] - 6.e7 * Y[1];
jac(1, 2) = -1.e4 * Y[1];
jac(2, 0) = 0.0;
jac(2, 1) = 6.e7 * Y[1];
jac(2, 2) = 0.0;
return jac;
}
std::vector<double>
backward_euler(const double t_start, const double t_end, const double dt_init,
const std::vector<double> y_init) {
const double tol{1.e-6};
const int max_iter{100};
int neq = static_cast<int>(y_init.size());
double time{t_start};
double dt{dt_init};
// starting point of the integration for each step
std::vector<double> y_n{y_init};
std::vector<double> y_new{0.0};
while (time < t_end) {
for (int n = 0; n < static_cast<int>(y_new.size()); ++n) {
y_new[n] = y_n[n];
}
bool converged{false};
int niter{0};
while (!converged && niter < max_iter) {
// construct the matrix A = I - dt J
Array A = jacobian(y_n);
for (int irow = 0; irow < neq; ++irow) {
for (int jcol = 0; jcol < neq; ++jcol) {
A(irow, jcol) *= -dt;
if (irow == jcol) {
A(irow, jcol) += 1.0;
}
}
}
// construct the HS
std::vector<double> b = rhs(y_new);
for (int irow = 0; irow < neq; ++irow) {
b[irow] = y_n[irow] - y_new[irow] + dt * b[irow];
}
// solve the linear system A dy = b
std::vector<double> dy = gauss_elim(A, b);
// correct our guess
for (int irow = 0; irow < neq; ++irow) {
y_new[irow] += dy[irow];
}
// check for convergence by computing the norms
double dy_norm{0.0};
double y_new_norm(0.0);
for (int irow = 0; irow < neq; ++irow) {
dy_norm += dy[irow] * dy[irow];
y_new_norm += y_new[irow] * y_new[irow];
}
double error = std::sqrt(dy_norm / y_new_norm);
if (error < tol) {
converged = true;
}
niter++;
}
if (!converged) {
std::cout << "convergence failure" << std::endl;
std::exit(EXIT_FAILURE);
}
if (time + dt > t_end) {
dt = t_end - time;
}
// reset the solution for the next step
for (int irow = 0; irow < neq; ++irow) {
y_n[irow] = y_new[irow];
}
time += dt;
}
return y_n;
}
int main() {
// setup the initial conditions
std::vector<double> y_init{1.0, 0.0, 0.0};
// we are going to do a bunch of short integrations to different
// end points -- this holds the end times for each interval
std::vector<double> t_ends{-6.0, -5.0, -4.0, -3.0, -2.0, -1.0, 0.0,
1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
double time{0.0};
std::vector<double> y_old{y_init};
for (auto log_t_max: t_ends) {
double tmax = std::pow(10.0, log_t_max);
auto y_new = backward_euler(time, tmax, tmax/10.0, y_old);
time = tmax;
std::cout << tmax << " " << y_new[0] << " " << y_new[1] << " " << y_new[2] << std::endl;
for (int irow = 0; irow < 3; ++irow) {
y_old[irow] = y_new[irow];
}
}
}
This can be compiled as:
g++ -std=c++17 -I. -o stiff_system stiff_system.cpp