13#ifndef SANDALS_RUNGEKUTTA_HH
14#define SANDALS_RUNGEKUTTA_HH
50 template <
typename Real, Integer S, Integer N, Integer M>
53 using VectorX = Eigen::Vector<Real, Eigen::Dynamic>;
54 using MatrixX = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
56 using MatrixK = Eigen::Matrix<Real, N, S>;
58 using MatrixJ = Eigen::Matrix<Real, N*S, N*S>;
59 using VectorP = Eigen::Matrix<Real, N+M, 1>;
60 using MatrixP = Eigen::Matrix<Real, N+M, N+M>;
61 using NewtonX = Optimist::RootFinder::Newton<Real, N>;
62 using NewtonK = Optimist::RootFinder::Newton<Real, N*S>;
77 using Time = Eigen::Vector<Real, Eigen::Dynamic>;
82 mutable Eigen::FullPivLU<MatrixP>
m_lu;
142 bool is_erk()
const {
return this->m_tableau.
type == Type::ERK;}
148 bool is_irk()
const {
return this->m_tableau.
type == Type::IRK;}
178 std::string
name()
const {
return this->m_tableau.
name;}
245 this->m_system = std::make_shared<ImplicitWrapper<Real, N, M>>(F, JF_x, JF_x_dot, h, Jh_x, in_domain);
267 this->m_system = std::make_shared<ImplicitWrapper<Real, N, M>>(
name, F, JF_x, JF_x_dot, h, Jh_x, in_domain);
285 this->m_system = std::make_shared<ExplicitWrapper<Real, N, M>>(f, Jf_x, h, Jh_x, in_domain);
305 this->m_system = std::make_shared<ExplicitWrapper<Real, N, M>>(
name, f, Jf_x, h, Jh_x, in_domain);
325 this->m_system = std::make_shared<LinearWrapper<Real, N, M>>(E,
A,
b, h, Jh_x, in_domain);
347 this->m_system = std::make_shared<LinearWrapper<Real, N, M>>(
name, E,
A,
b, h, Jh_x, in_domain);
369 this->m_system = std::make_shared<SemiExplicitWrapper<Real, N, M>>(
A, TA_x,
b, Jb_x, h, Jh_x, in_domain);
393 this->m_system = std::make_shared<SemiExplicitWrapper<Real, N, M>>(
name,
A, TA_x,
b, Jb_x, h, Jh_x, in_domain);
412 void absolute_tolerance(Real
const t_absolute_tolerance) {this->m_absolute_tolerance = t_absolute_tolerance;}
424 void relative_tolerance(Real
const t_relative_tolerance) {this->m_relative_tolerance = t_relative_tolerance;}
436 void safety_factor(Real
const t_safety_factor) {this->m_safety_factor = t_safety_factor;}
448 void min_safety_factor(Real
const t_min_safety_factor) {this->m_min_safety_factor = t_min_safety_factor;}
460 void max_safety_factor(Real
const t_max_safety_factor) {this->m_max_safety_factor = t_max_safety_factor;}
472 void min_step(Real
const t_min_step) {this->m_min_step = t_min_step;}
496 void adaptive(
bool t_adaptive) {this->m_adaptive = t_adaptive;}
519 this->m_verbose = t_verbose;
520 this->m_newtonX.verbose_mode(t_verbose);
521 this->m_newtonK.verbose_mode(t_verbose);
544 void reverse(
bool t_reverse) {this->m_reverse = t_reverse;}
579 {this->m_projection_tolerance = t_projection_tolerance;}
592 {this->m_max_projection_iterations = t_max_projection_iterations;}
604 void projection(
bool t_projection) {this->m_projection = t_projection;}
652 Real desired_error{this->m_absolute_tolerance + this->m_relative_tolerance *
653 std::max(x.array().abs().maxCoeff(), x_e.array().abs().maxCoeff())};
654 Real truncation_error{(x - x_e).array().abs().maxCoeff()};
655 return h_k * std::min(this->m_max_safety_factor, std::max(this->m_min_safety_factor,
656 this->m_safety_factor * std::pow(desired_error/truncation_error,
657 1.0/std::max(this->m_tableau.
order, this->m_tableau.order_e))));
666 std::ostringstream os;
668 <<
"Runge-Kutta method:\t" << this->
name() << std::endl
669 <<
"\t- order:\t" << this->
order() << std::endl
670 <<
"\t- stages:\t" << this->
stages() << std::endl
672 switch (this->
type()) {
673 case Type::ERK: os <<
"explicit";
break;
674 case Type::IRK: os <<
"implicit";
break;
675 case Type::DIRK: os <<
"diagonally implicit";
break;
679 <<
"\t- embedded:\t" << this->
is_embedded() << std::endl;
681 os <<
"\t- system:\t" << this->m_system->name() << std::endl;
683 os <<
"\t- system:\t" <<
"none" << std::endl;
718 Real & h_new,
MatrixK & K)
const
725 for (
Integer i{0}; i < S; ++i) {
726 x_node = x_old + K(all, seqN(0, i)) * this->m_tableau.
A(i, seqN(0, i)).transpose();
727 if (!this->m_reverse) {
728 K.col(i) = h_old *
static_cast<Explicit<Real, N, M> const *
>(this->m_system.get())->
f(x_node, t_old + h_old*this->m_tableau.
c(i));
733 if (!K.allFinite()) {
return false;}
736 x_new = x_old + K * this->m_tableau.
b;
739 if (this->m_adaptive && this->m_tableau.
is_embedded) {
740 VectorN x_emb = x_old + K * this->m_tableau.
b_e;
771 VectorN x_node(x + K(all, seqN(0, s)) * this->m_tableau.
A(s, seqN(0, s)).transpose());
772 if (!this->m_reverse) {
773 fun = this->m_system->F(x_node, K.col(s)/h, t + h * this->m_tableau.c(s));
775 fun = this->m_system->F_reverse(x_node, K.col(s)/h, t + h * this->m_tableau.c(s));
815 VectorN x_node(x + K(all, seqN(0, s)) * this->m_tableau.
A(s, seqN(0, s)).transpose());
816 if (!this->m_reverse) {
817 jac = this->m_system->JF_x_dot(x_node, K.col(s)/h, t + h * this->m_tableau.c(s)) / h;
819 jac = this->m_system->JF_x_dot_reverse(x_node, K.col(s)/h, t + h * this->m_tableau.c(s)) / h;
838 Real & h_new,
MatrixK & K)
const
841 VectorN K_ini(VectorN::Zero());
844 for (
Integer s{0}; s < S; ++s) {
845 if (this->m_newtonX.solve(
846 [
this, s, &K, &x_old, t_old, h_old](
VectorN const & K_fun,
VectorN & fun)
847 {K.col(s) = K_fun; this->erk_implicit_function(s, x_old, t_old, h_old, K, fun);},
848 [
this, s, &K, &x_old, t_old, h_old](
VectorN const & K_jac,
MatrixN & jac)
849 {K.col(s) = K_jac; this->erk_implicit_jacobian(s, x_old, t_old, h_old, K, jac);},
858 x_new = x_old + K * this->m_tableau.
b;
861 if (this->m_adaptive && this->m_tableau.
is_embedded) {
862 VectorN x_emb(x_old + K * this->m_tableau.
b_e);
931 MatrixK K_mat{K.reshaped(N, S)};
933 for (
Integer i{0}; i < S; ++i) {
934 x_node = x + K_mat * this->m_tableau.
A.row(i).transpose();
935 if (!this->m_reverse) {
936 fun_mat.col(i) = this->m_system->F(x_node, K_mat.col(i)/h, t + h * this->m_tableau.c(i));
938 fun_mat.col(i) = this->m_system->F_reverse(x_node, K_mat.col(i)/h, t + h * this->m_tableau.c(i));
941 fun = fun_mat.reshaped(N*S, 1);
988 MatrixK K_mat{K.reshaped(N, S)};
992 auto idx = seqN(0, N), jdx = seqN(0, N);
993 for (
Integer i{0}; i < S; ++i) {
994 t_node = t + h * this->m_tableau.
c(i);
995 x_node = x + K_mat * this->m_tableau.
A.row(i).transpose();
998 x_dot_node = K_mat.col(i) / h;
999 if (!this->m_reverse) {
1000 JF_x = this->m_system->JF_x(x_node, x_dot_node, t_node);
1001 JF_x_dot = this->m_system->JF_x_dot(x_node, x_dot_node, t_node);
1003 JF_x = this->m_system->JF_x_reverse(x_node, x_dot_node, t_node);
1004 JF_x_dot = this->m_system->JF_x_dot_reverse(x_node, x_dot_node, t_node);
1009 for (
Integer j{0}; j < S; ++j) {
1012 jac(idx, jdx) = this->m_tableau.
A(i,j) * JF_x + JF_x_dot / h;
1014 jac(idx, jdx) = this->m_tableau.
A(i,j) * JF_x;
1035 Real & h_new,
MatrixK & K)
const
1038 VectorK K_ini(VectorK::Zero());
1041 if (!this->m_newtonK.solve(
1042 [
this, &x_old, t_old, h_old](
VectorK const & K_fun,
VectorK & fun)
1043 {this->irk_function(x_old, t_old, h_old, K_fun, fun);},
1044 [
this, &x_old, t_old, h_old](
VectorK const & K_jac,
MatrixJ & jac)
1045 {this->irk_jacobian(x_old, t_old, h_old, K_jac, jac);},
1050 K = K_vec.reshaped(N, S);
1051 if (!K.allFinite()) {
return false;}
1054 x_new = x_old + K * this->m_tableau.
b;
1057 if (this->m_adaptive && this->m_tableau.
is_embedded) {
1058 VectorN x_emb(x_old + K * this->m_tableau.
b_e);
1098 VectorN x_node(x + K(all, seqN(0, n+1)) * this->m_tableau.
A(n, seqN(0, n+1)).transpose());
1099 if (!this->m_reverse) {
1100 fun = this->m_system->F(x_node, K.col(n)/h, t + h * this->m_tableau.c(n));
1102 fun = this->m_system->F_reverse(x_node, K.col(n)/h, t + h * this->m_tableau.c(n));
1144 Real t_node{t + h * this->m_tableau.
c(n)};
1145 VectorN x_node(x + K(all, seqN(0, n+1)) * this->m_tableau.
A(n, seqN(0, n+1)).transpose());
1146 VectorN x_dot_node(K.col(n)/h);
1147 if (!this->m_reverse) {
1148 jac = this->m_tableau.
A(n,n) * this->m_system->JF_x(x_node, x_dot_node, t_node) +
1149 this->m_system->JF_x_dot(x_node, x_dot_node, t_node) / h;
1151 jac = this->m_tableau.
A(n,n) * this->m_system->JF_x_reverse(x_node, x_dot_node, t_node) +
1152 this->m_system->JF_x_dot_reverse(x_node, x_dot_node, t_node) / h;
1171 Real & h_new,
MatrixK & K)
const
1174 VectorN K_ini(VectorN::Zero());
1177 for (
Integer n{0}; n < S; ++n) {
1178 if (this->m_newtonX.solve(
1179 [
this, n, &K, &x_old, t_old, h_old](
VectorN const & K_fun,
VectorN & fun)
1180 {K.col(n) = K_fun; this->dirk_function(n, x_old, t_old, h_old, K, fun);},
1181 [
this, n, &K, &x_old, t_old, h_old](
VectorN const & K_jac,
MatrixN & jac)
1182 {K.col(n) = K_jac; this->dirk_jacobian(n, x_old, t_old, h_old, K, jac);},
1191 x_new = x_old + K * this->m_tableau.
b;
1194 if (this->m_adaptive && this->m_tableau.
is_embedded) {
1195 VectorN x_emb(x_old + K * this->m_tableau.
b_e);
1216 #define CMD "Sandals::RungeKutta::step(...): "
1218 SANDALS_ASSERT(this->m_system->in_domain(x_old, t_old),
CMD "in " << this->m_tableau.name <<
1219 " solver, at t = " << t_old <<
", x = " << x_old.transpose() <<
", system out of domain.");
1221 if (this->
is_erk() && this->m_system->is_explicit()) {
1223 }
else if (this->
is_erk() && this->m_system->is_implicit()) {
1226 return this->
dirk_step(x_old, t_old, h_old, x_new, h_new, K);
1228 return this->
irk_step(x_old, t_old, h_old, x_new, h_new, K);
1249 #define CMD "Sandals::RungeKutta::advance(...): "
1253 h_old <<
", expected > 0.");
1257 if (!this->
step(x_old, t_old, h_old, x_new, h_new, K))
1260 Real t_tmp{t_old}, h_tmp{h_old / Real(2.0)};
1263 Integer max_k{this->m_max_substeps * this->m_max_substeps}, k{2};
1267 if (this->
step(x_tmp, t_tmp, h_tmp, x_new, h_new_tmp, K)) {
1273 if (k > 0 && k < max_k) {
1277 h_tmp = Real(2.0) * h_tmp;
1278 if (this->m_verbose) {
1280 ", integration succedded disable one substepping layer.");
1286 SANDALS_ASSERT(std::isfinite(x_tmp.maxCoeff()),
CMD "in " << this->m_tableau.name <<
1287 " solver, at t = " << t_tmp <<
", ||x||_inf = inf, computation interrupted.");
1294 t_tmp <<
", integration failed with h = " << h_tmp <<
", aborting.");
1299 "at t = " << t_tmp <<
", integration failed, adding substepping layer.");}
1314 if (this->m_projection) {
1316 if (this->
project(x_new, t_old + h_new, x_projected)) {
1317 x_new = x_projected;
1343 sol.
resize(t_mesh.size());
1346 sol.
t(0) = t_mesh(0);
1348 sol.
h.col(0) = this->m_system->h(ics, t_mesh(0));
1351 if (this->m_step_callback) {this->
m_step_callback(0, ics, t_mesh(0));}
1355 VectorN x_old_step(ics), x_new_step(ics);
1356 Real t_step{t_mesh(0)}, h_step{t_mesh(1) - t_mesh(0)}, h_tmp_step{h_step}, h_new_step;
1357 bool mesh_point_bool, saturation_bool;
1361 if (!this->
advance(x_old_step, t_step, h_step, x_new_step, h_new_step)) {
return false;}
1369 if (this->m_adaptive && this->m_tableau.
is_embedded && !mesh_point_bool && saturation_bool) {
1370 h_tmp_step = h_new_step;
1371 h_step = t_mesh(
step+1) - t_step;
1373 h_step = h_new_step;
1377 if (!this->m_adaptive || mesh_point_bool) {
1381 h_step = h_tmp_step;
1384 sol.
t(
step) = t_step;
1385 sol.
x.col(
step) = x_new_step;
1386 sol.
h.col(
step) = this->m_system->h(x_new_step, t_step);
1389 if (this->m_step_callback) {this->
m_step_callback(step, x_new_step, t_step);}
1392 if (std::abs(t_step - t_mesh(last)) <
SQRT_EPSILON) {
break;}
1395 x_old_step = x_new_step;
1417 #define CMD "Sandals::RungeKutta::adaptive_solve(...): "
1422 return this->
solve(t_mesh, ics, sol);
1423 }
else if (!this->m_adaptive) {
1425 return this->
solve(t_mesh, ics, sol);
1429 Real t_step{t_mesh(0)}, h_step{t_mesh(1) - t_mesh(0)}, h_new_step, scale{100.0};
1430 Real h_min{std::max(this->m_min_step, h_step/scale)}, h_max{scale*h_step};
1432 Integer safety_length{
static_cast<Integer>(std::ceil(std::abs(t_mesh(last) - t_mesh(0))/(2.0*h_min)))};
1433 sol.
resize(safety_length);
1435 sol.
resize(t_mesh.size());
1439 sol.
t(0) = t_mesh(0);
1441 sol.
h.col(0) = this->m_system->h(ics, t_mesh(0));
1444 if (this->m_step_callback) {this->
m_step_callback(0, ics, t_mesh(0));}
1448 VectorN x_old_step(ics), x_new_step(ics);
1452 this->
advance(x_old_step, t_step, h_step, x_new_step, h_new_step);
1458 if (this->m_adaptive && this->m_tableau.
is_embedded) {
1459 h_step = std::max(std::min(h_new_step, h_max), h_min);
1468 sol.
t(
step) = t_step;
1469 sol.
x.col(
step) = x_new_step;
1470 sol.
h.col(
step) = this->m_system->h(x_new_step, t_step);
1473 if (this->m_step_callback) {this->
m_step_callback(step, x_new_step, t_step);}
1476 if (std::abs(t_step - t_mesh(last)) <
SQRT_EPSILON) {
break;}
1477 else if (t_step + h_step > t_mesh(last)) {h_step = t_mesh(last) - t_step;}
1480 x_old_step = x_new_step;
1502 #define CMD "Sandals::RungeKutta::project(...): "
1513 A.template block<N, N>(0, 0) = MatrixN::Identity();
1514 for (
Integer k{0}; k < this->m_max_projection_iterations; ++k) {
1524 h = this->m_system->h(x_projected, t);
1525 Jh_x = this->m_system->Jh_x(x_projected, t);
1528 if (h.norm() < this->m_projection_tolerance) {
return true;}
1531 A.template block<N, M>(0, N) = Jh_x.transpose();
1532 A.template block<M, N>(N, 0) = Jh_x;
1533 b.template head<N>() = x - x_projected;
1534 b.template tail<M>() = -h;
1537 this->m_lu.compute(
A);
1539 x_step = this->m_lu.solve(
b);
1542 if (x_step.norm() < this->m_projection_tolerance * this->m_projection_tolerance) {
return false;}
1545 x_projected.noalias() += x_step(Eigen::seqN(0, N));
1568 std::vector<Integer>
const & projected_invariants,
VectorN & x_projected)
const
1570 #define CMD "Sandals::RungeKutta::project_ics(...): "
1584 A.block(0, 0, X, X) = MatrixX::Identity(X+H, X+H);
1585 Eigen::FullPivLU<MatrixX> lu;
1586 for (
Integer k{0}; k < this->m_max_projection_iterations; ++k) {
1596 h = this->m_system->h(x_projected, t);
1597 Jh_x = this->m_system->Jh_x(x_projected, t);
1600 h = h(projected_invariants);
1601 Jh_x = Jh_x(projected_invariants, projected_equations);
1604 if (h.norm() < this->m_projection_tolerance) {
return true;}
1607 A.block(0, X, X, H) = Jh_x.transpose();
1608 A.block(X, 0, H, X) = Jh_x;
1609 b.head(X) = x(projected_equations) - x_projected(projected_equations);
1615 x_step = this->m_lu.solve(
b);
1618 if (x_step.norm() < this->m_projection_tolerance * this->m_projection_tolerance) {
return false;}
1621 x_projected(projected_equations).noalias() += x_step;
1643 #define CMD "Sandals::RungeKutta::estimate_order(...): "
1647 for (
Integer i{0}; i < static_cast<Integer>(t_mesh.size()); ++i) {
1651 CMD "expected the same initial time.");
1653 CMD "expected the same final time.");
1656 for (
Integer j{1}; j < static_cast<Integer>(t_mesh[i].size()); ++j) {
1658 CMD "expected a fixed step.");
1665 VectorX h_vec(t_mesh.size()), e_vec(t_mesh.size());
1666 for (
Integer i{0}; i < static_cast<Integer>(t_mesh.size()); ++i) {
1668 "the" << i <<
"-th time mesh.");
1669 sol_ana = sol(sol_num.
t);
1671 CMD "expected the same number of states in analytical solution.");
1673 CMD "expected the same number of steps in analytical solution.");
1674 h_vec(i) = std::abs(sol_num.
t(1) - sol_num.
t(0));
1675 e_vec(i) = (sol_ana - sol_num.
x).array().abs().maxCoeff();
1681 return ((
A.transpose() *
A).ldlt().solve(
A.transpose() *
b))(0);
#define SANDALS_BASIC_CONSTANTS(Real)
Definition Sandals.hh:70
#define SANDALS_ASSERT(COND, MSG)
Definition Sandals.hh:44
#define SANDALS_WARNING(MSG)
Definition Sandals.hh:53
Class container for the system of explicit ODEs.
Definition Explicit.hh:42
virtual VectorF f(VectorF const &x, Real const t) const =0
VectorF f_reverse(VectorF const &x, Real const t) const
Definition Explicit.hh:158
static const FunctionH DefaultH
Definition Explicit.hh:254
std::function< MatrixJF(VectorF const &, Real const)> FunctionJF
Definition Explicit.hh:249
static const FunctionID DefaultID
Definition Explicit.hh:256
std::function< VectorH(VectorF const &, Real const)> FunctionH
Definition Explicit.hh:250
std::function< bool(VectorF const &, Real const)> FunctionID
Definition Explicit.hh:252
static const FunctionJH DefaultJH
Definition Explicit.hh:255
std::function< MatrixJH(VectorF const &, Real const)> FunctionJH
Definition Explicit.hh:251
std::function< VectorF(VectorF const &, Real const)> FunctionF
Definition Explicit.hh:248
Eigen::Matrix< Real, N, N > MatrixJF
Definition Implicit.hh:47
Eigen::Vector< Real, N > VectorF
Definition Implicit.hh:46
Eigen::Vector< Real, M > VectorH
Definition Implicit.hh:48
std::shared_ptr< Implicit< Real, N, M > > Pointer
Definition Implicit.hh:45
Eigen::Matrix< Real, M, N > MatrixJH
Definition Implicit.hh:49
std::function< bool(VectorF const &, Real const)> FunctionID
Definition Implicit.hh:279
std::function< VectorF(VectorF const &, VectorF const &, Real const)> FunctionF
Definition Implicit.hh:275
std::function< MatrixJF(VectorF const &, VectorF const &, Real const)> FunctionJF
Definition Implicit.hh:276
std::function< VectorH(VectorF const &, Real const)> FunctionH
Definition Implicit.hh:277
std::function< MatrixJH(VectorF const &, Real const)> FunctionJH
Definition Implicit.hh:278
static const FunctionID DefaultID
Definition Implicit.hh:283
static const FunctionH DefaultH
Definition Implicit.hh:281
static const FunctionJH DefaultJH
Definition Implicit.hh:282
std::function< VectorB(Real const)> FunctionB
Definition Linear.hh:212
static const FunctionH DefaultH
Definition Linear.hh:217
std::function< MatrixE(Real const)> FunctionE
Definition Linear.hh:210
static const FunctionJH DefaultJH
Definition Linear.hh:218
static const FunctionID DefaultID
Definition Linear.hh:219
std::function< MatrixJH(VectorF const &, Real const)> FunctionJH
Definition Linear.hh:214
std::function< bool(VectorF const &, Real const)> FunctionID
Definition Linear.hh:215
std::function< MatrixA(Real const)> FunctionA
Definition Linear.hh:211
std::function< VectorH(VectorF const &, Real const)> FunctionH
Definition Linear.hh:213
Real & min_step()
Definition RungeKutta.hh:466
bool projection()
Definition RungeKutta.hh:598
bool adaptive_mode()
Definition RungeKutta.hh:490
bool is_erk() const
Definition RungeKutta.hh:142
typename Tableau< Real, S >::Matrix MatrixS
Definition RungeKutta.hh:64
Real estimate_order(std::vector< VectorX > const &t_mesh, VectorN const &ics, std::function< MatrixX(VectorX)> &sol) const
Definition RungeKutta.hh:1639
void enable_reverse_mode()
Definition RungeKutta.hh:549
void projection(bool t_projection)
Definition RungeKutta.hh:604
void dirk_function(Integer n, VectorN const &x, Real const t, Real const h, MatrixK const &K, VectorN &fun) const
Definition RungeKutta.hh:1093
typename Tableau< Real, S >::Vector VectorS
Definition RungeKutta.hh:63
Real & min_safety_factor()
Definition RungeKutta.hh:442
void semi_explicit_system(typename SemiExplicitWrapper< Real, N, M >::FunctionA A, typename SemiExplicitWrapper< Real, N, M >::FunctionTA TA_x, typename SemiExplicitWrapper< Real, N, M >::FunctionB b, typename SemiExplicitWrapper< Real, N, M >::FunctionJB Jb_x, typename SemiExplicitWrapper< Real, N, M >::FunctionH h=SemiExplicitWrapper< Real, N, M >::DefaultH, typename SemiExplicitWrapper< Real, N, M >::FunctionJH Jh_x=SemiExplicitWrapper< Real, N, M >::DefaultJH, typename SemiExplicitWrapper< Real, N, M >::FunctionID in_domain=SemiExplicitWrapper< Real, N, M >::DefaultID)
Definition RungeKutta.hh:360
bool reverse_mode()
Definition RungeKutta.hh:538
typename Implicit< Real, N, M >::MatrixJF MatrixN
Definition RungeKutta.hh:66
void erk_implicit_function(Integer const s, VectorN const &x, Real const t, Real const h, MatrixK const &K, VectorN &fun) const
Definition RungeKutta.hh:766
Real & safety_factor()
Definition RungeKutta.hh:430
Real m_absolute_tolerance
Definition RungeKutta.hh:86
Real m_min_safety_factor
Definition RungeKutta.hh:89
RungeKutta & operator=(RungeKutta const &)=delete
typename Implicit< Real, N, M >::MatrixJH MatrixM
Definition RungeKutta.hh:68
Eigen::Matrix< Real, N+M, 1 > VectorP
Definition RungeKutta.hh:59
Tableau< Real, S > & tableau()
Definition RungeKutta.hh:160
System system()
Definition RungeKutta.hh:220
VectorS b_embedded() const
Definition RungeKutta.hh:208
RungeKutta(const RungeKutta &)=delete
FunctionSC step_callback()
Definition RungeKutta.hh:560
bool advance(VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new) const
Definition RungeKutta.hh:1247
typename Implicit< Real, N, M >::Pointer System
Definition RungeKutta.hh:75
bool step(VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new, MatrixK &K) const
Definition RungeKutta.hh:1213
void enable_projection()
Definition RungeKutta.hh:609
Eigen::Matrix< Real, N *N, S > MatrixJK
Definition RungeKutta.hh:57
RungeKutta(Tableau< Real, S > const &t_tableau, System t_system)
Definition RungeKutta.hh:127
Real m_relative_tolerance
Definition RungeKutta.hh:87
Integer & max_substeps()
Definition RungeKutta.hh:478
typename Implicit< Real, N, M >::VectorF VectorN
Definition RungeKutta.hh:65
void linear_system(std::string name, typename LinearWrapper< Real, N, M >::FunctionE E, typename LinearWrapper< Real, N, M >::FunctionA A, typename LinearWrapper< Real, N, M >::FunctionB b, typename LinearWrapper< Real, N, M >::FunctionH h=LinearWrapper< Real, N, M >::DefaultH, typename LinearWrapper< Real, N, M >::FunctionJH Jh_x=LinearWrapper< Real, N, M >::DefaultJH, typename LinearWrapper< Real, N, M >::FunctionID in_domain=LinearWrapper< Real, N, M >::DefaultID)
Definition RungeKutta.hh:338
RungeKutta(Tableau< Real, S > const &t_tableau)
Definition RungeKutta.hh:117
Eigen::Matrix< Real, N+M, N+M > MatrixP
Definition RungeKutta.hh:60
void adaptive(bool t_adaptive)
Definition RungeKutta.hh:496
Integer order() const
Definition RungeKutta.hh:184
void projection_tolerance(Real const t_projection_tolerance)
Definition RungeKutta.hh:578
Eigen::Matrix< Real, N *S, N *S > MatrixJ
Definition RungeKutta.hh:58
void irk_jacobian(VectorN const &x, Real const t, Real const h, VectorK const &K, MatrixJ &jac) const
Definition RungeKutta.hh:980
VectorS b() const
Definition RungeKutta.hh:202
typename Tableau< Real, S >::Type Type
Definition RungeKutta.hh:76
bool erk_explicit_step(VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new, MatrixK &K) const
Definition RungeKutta.hh:717
NewtonX m_newtonX
Definition RungeKutta.hh:80
bool is_dirk() const
Definition RungeKutta.hh:154
VectorS c() const
Definition RungeKutta.hh:214
void reverse(bool t_reverse)
Definition RungeKutta.hh:544
bool m_reverse
Definition RungeKutta.hh:95
bool m_adaptive
Definition RungeKutta.hh:93
void max_safety_factor(Real const t_max_safety_factor)
Definition RungeKutta.hh:460
Eigen::Vector< Real, Eigen::Dynamic > VectorX
Definition RungeKutta.hh:53
std::string name() const
Definition RungeKutta.hh:178
Real m_min_step
Definition RungeKutta.hh:91
void implicit_system(typename ImplicitWrapper< Real, N, M >::FunctionF F, typename ImplicitWrapper< Real, N, M >::FunctionJF JF_x, typename ImplicitWrapper< Real, N, M >::FunctionJF JF_x_dot, typename ImplicitWrapper< Real, N, M >::FunctionH h=ImplicitWrapper< Real, N, M >::DefaultH, typename ImplicitWrapper< Real, N, M >::FunctionJH Jh_x=ImplicitWrapper< Real, N, M >::DefaultJH, typename ImplicitWrapper< Real, N, M >::FunctionID in_domain=ImplicitWrapper< Real, N, M >::DefaultID)
Definition RungeKutta.hh:237
void semi_explicit_system(std::string name, typename SemiExplicitWrapper< Real, N, M >::FunctionA A, typename SemiExplicitWrapper< Real, N, M >::FunctionTA TA_x, typename SemiExplicitWrapper< Real, N, M >::FunctionB b, typename SemiExplicitWrapper< Real, N, M >::FunctionJB Jb_x, typename SemiExplicitWrapper< Real, N, M >::FunctionH h=SemiExplicitWrapper< Real, N, M >::DefaultH, typename SemiExplicitWrapper< Real, N, M >::FunctionJH Jh_x=SemiExplicitWrapper< Real, N, M >::DefaultJH, typename SemiExplicitWrapper< Real, N, M >::FunctionID in_domain=SemiExplicitWrapper< Real, N, M >::DefaultID)
Definition RungeKutta.hh:383
void enable_adaptive_mode()
Definition RungeKutta.hh:501
Real projection_tolerance()
Definition RungeKutta.hh:572
void verbose_mode(bool t_verbose)
Definition RungeKutta.hh:518
std::function< void(Integer const, VectorX const &, Real const)> FunctionSC
Definition RungeKutta.hh:69
void absolute_tolerance(Real const t_absolute_tolerance)
Definition RungeKutta.hh:412
void erk_implicit_jacobian(Integer const s, VectorN const &x, Real const t, Real const h, MatrixK const &K, MatrixN &jac) const
Definition RungeKutta.hh:810
Real relative_tolerance()
Definition RungeKutta.hh:418
bool verbose_mode()
Definition RungeKutta.hh:512
void linear_system(typename LinearWrapper< Real, N, M >::FunctionE E, typename LinearWrapper< Real, N, M >::FunctionA A, typename LinearWrapper< Real, N, M >::FunctionB b, typename LinearWrapper< Real, N, M >::FunctionH h=LinearWrapper< Real, N, M >::DefaultH, typename LinearWrapper< Real, N, M >::FunctionJH Jh_x=LinearWrapper< Real, N, M >::DefaultJH, typename LinearWrapper< Real, N, M >::FunctionID in_domain=LinearWrapper< Real, N, M >::DefaultID)
Definition RungeKutta.hh:317
const Real SQRT_EPSILON
Definition RungeKutta.hh:73
void explicit_system(typename ExplicitWrapper< Real, N, M >::FunctionF f, typename ExplicitWrapper< Real, N, M >::FunctionJF Jf_x, typename ExplicitWrapper< Real, N, M >::FunctionH h=ExplicitWrapper< Real, N, M >::DefaultH, typename ExplicitWrapper< Real, N, M >::FunctionJH Jh_x=ExplicitWrapper< Real, N, M >::DefaultJH, typename ExplicitWrapper< Real, N, M >::FunctionID in_domain=ExplicitWrapper< Real, N, M >::DefaultID)
Definition RungeKutta.hh:278
Eigen::FullPivLU< MatrixP > m_lu
Definition RungeKutta.hh:82
System m_system
Definition RungeKutta.hh:85
void min_safety_factor(Real const t_min_safety_factor)
Definition RungeKutta.hh:448
Tableau< Real, S > m_tableau
Definition RungeKutta.hh:84
bool solve(VectorX const &t_mesh, VectorN const &ics, Solution< Real, N, M > &sol) const
Definition RungeKutta.hh:1338
bool adaptive_solve(VectorX const &t_mesh, VectorN const &ics, Solution< Real, N, M > &sol) const
Definition RungeKutta.hh:1412
void relative_tolerance(Real const t_relative_tolerance)
Definition RungeKutta.hh:424
void max_substeps(Integer const t_max_substeps)
Definition RungeKutta.hh:484
void min_step(Real const t_min_step)
Definition RungeKutta.hh:472
Eigen::Matrix< Real, N, S > MatrixK
Definition RungeKutta.hh:56
Integer stages() const
Definition RungeKutta.hh:172
void enable_verbose_mode()
Definition RungeKutta.hh:527
Integer & max_projection_iterations()
Definition RungeKutta.hh:585
bool is_irk() const
Definition RungeKutta.hh:148
Optimist::RootFinder::Newton< Real, N > NewtonX
Definition RungeKutta.hh:61
void step_callback(FunctionSC const &t_step_callback)
Definition RungeKutta.hh:566
bool erk_implicit_step(VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new, MatrixK &K) const
Definition RungeKutta.hh:837
Real m_projection_tolerance
Definition RungeKutta.hh:98
FunctionSC m_step_callback
Definition RungeKutta.hh:96
Real estimate_step(VectorN const &x, VectorN const &x_e, Real const h_k) const
Definition RungeKutta.hh:650
void dirk_jacobian(Integer n, VectorN const &x, Real const t, Real const h, MatrixK const &K, MatrixN &jac) const
Definition RungeKutta.hh:1139
void system(System t_system)
Definition RungeKutta.hh:226
void disable_adaptive_mode()
Definition RungeKutta.hh:506
Tableau< Real, S > const & tableau() const
Definition RungeKutta.hh:166
Real m_max_safety_factor
Definition RungeKutta.hh:90
bool m_verbose
Definition RungeKutta.hh:94
Real m_safety_factor
Definition RungeKutta.hh:88
bool project_ics(VectorN const &x, Real const t, std::vector< Integer > const &projected_equations, std::vector< Integer > const &projected_invariants, VectorN &x_projected) const
Definition RungeKutta.hh:1567
void disable_reverse_mode()
Definition RungeKutta.hh:554
bool dirk_step(VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new, MatrixK &K) const
Definition RungeKutta.hh:1170
Eigen::Vector< Real, N *S > VectorK
Definition RungeKutta.hh:55
void info(std::ostream &os)
Definition RungeKutta.hh:692
Type type() const
Definition RungeKutta.hh:136
void disable_projection()
Definition RungeKutta.hh:614
Eigen::Matrix< Real, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Definition RungeKutta.hh:54
void max_projection_iterations(Integer const t_max_projection_iterations)
Definition RungeKutta.hh:591
void safety_factor(Real const t_safety_factor)
Definition RungeKutta.hh:436
Eigen::Vector< Real, Eigen::Dynamic > Time
Definition RungeKutta.hh:77
void disable_verbose_mode()
Definition RungeKutta.hh:532
Integer m_max_substeps
Definition RungeKutta.hh:92
bool is_embedded() const
Definition RungeKutta.hh:190
std::string info() const
Definition RungeKutta.hh:664
void explicit_system(std::string name, typename ExplicitWrapper< Real, N, M >::FunctionF f, typename ExplicitWrapper< Real, N, M >::FunctionJF Jf_x, typename ExplicitWrapper< Real, N, M >::FunctionH h=ExplicitWrapper< Real, N, M >::DefaultH, typename ExplicitWrapper< Real, N, M >::FunctionJH Jh_x=ExplicitWrapper< Real, N, M >::DefaultJH, typename ExplicitWrapper< Real, N, M >::FunctionID in_domain=ExplicitWrapper< Real, N, M >::DefaultID)
Definition RungeKutta.hh:297
bool has_system()
Definition RungeKutta.hh:400
void irk_function(VectorN const &x, Real const t, Real const h, VectorK const &K, VectorK &fun) const
Definition RungeKutta.hh:928
typename Implicit< Real, N, M >::VectorH VectorM
Definition RungeKutta.hh:67
MatrixS A() const
Definition RungeKutta.hh:196
Real absolute_tolerance()
Definition RungeKutta.hh:406
Optimist::RootFinder::Newton< Real, N *S > NewtonK
Definition RungeKutta.hh:62
Integer m_max_projection_iterations
Definition RungeKutta.hh:99
bool project(VectorN const &x, Real const t, VectorN &x_projected) const
Definition RungeKutta.hh:1500
NewtonK m_newtonK
Definition RungeKutta.hh:81
void implicit_system(std::string name, typename ImplicitWrapper< Real, N, M >::FunctionF F, typename ImplicitWrapper< Real, N, M >::FunctionJF JF_x, typename ImplicitWrapper< Real, N, M >::FunctionJF JF_x_dot, typename ImplicitWrapper< Real, N, M >::FunctionH h=ImplicitWrapper< Real, N, M >::DefaultH, typename ImplicitWrapper< Real, N, M >::FunctionJH Jh_x=ImplicitWrapper< Real, N, M >::DefaultJH, typename ImplicitWrapper< Real, N, M >::FunctionID in_domain=ImplicitWrapper< Real, N, M >::DefaultID)
Definition RungeKutta.hh:258
Real & max_safety_factor()
Definition RungeKutta.hh:454
bool irk_step(VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new, MatrixK &K) const
Definition RungeKutta.hh:1034
bool m_projection
Definition RungeKutta.hh:100
std::function< MatrixJB(VectorF const &, Real const)> FunctionJB
Definition SemiExplicit.hh:273
static const FunctionH DefaultH
Definition SemiExplicit.hh:278
std::function< bool(VectorF const &, Real const)> FunctionID
Definition SemiExplicit.hh:276
std::function< VectorH(VectorF const &, Real const)> FunctionH
Definition SemiExplicit.hh:274
std::function< VectorB(VectorF const &, Real const)> FunctionB
Definition SemiExplicit.hh:272
std::function< TensorTA(VectorF const &, Real const)> FunctionTA
Definition SemiExplicit.hh:271
std::function< MatrixA(VectorF const &, Real const)> FunctionA
Definition SemiExplicit.hh:270
static const FunctionJH DefaultJH
Definition SemiExplicit.hh:279
static const FunctionID DefaultID
Definition SemiExplicit.hh:280
std::function< MatrixJH(VectorF const &, Real const)> FunctionJH
Definition SemiExplicit.hh:275
The namespace for the Sandals library.
Definition Sandals.hh:89
SANDALS_DEFAULT_INTEGER_TYPE Integer
The Integer type as used for the API.
Definition Sandals.hh:97
Class container for the numerical solution of a system of ODEs/DAEs.
Definition Solution.hh:62
MatrixN x
Definition Solution.hh:68
Integer size() const
Definition Solution.hh:126
void resize(Integer const size)
Definition Solution.hh:87
void conservative_resize(Integer const size)
Definition Solution.hh:97
Vector t
Definition Solution.hh:67
MatrixM h
Definition Solution.hh:69
Struct container for the Butcher tableau of a Runge-Kutta method.
Definition Tableau.hh:38
enum class type :Integer {ERK=0, IRK=1, DIRK=2} Type
Definition Tableau.hh:42
Type type
Definition Tableau.hh:47
Integer order
Definition Tableau.hh:48
std::string name
Definition Tableau.hh:46
Vector b_e
Definition Tableau.hh:52
Eigen::Matrix< Real, S, S > Matrix
Definition Tableau.hh:44
Matrix A
Definition Tableau.hh:50
Eigen::Vector< Real, S > Vector
Definition Tableau.hh:43
Vector c
Definition Tableau.hh:53
Vector b
Definition Tableau.hh:51
bool is_embedded
Definition Tableau.hh:54