13#ifndef SANDALS_RUNGEKUTTA_HH
14#define SANDALS_RUNGEKUTTA_HH
48 template <
typename Real, Integer S, Integer N, Integer M>
51 using VectorX = Eigen::Vector<Real, Eigen::Dynamic>;
52 using MatrixX = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
54 using MatrixK = Eigen::Matrix<Real, N, S>;
55 using MatrixJ = Eigen::Matrix<Real, N*S, N*S>;
56 using VectorP = Eigen::Matrix<Real, N+M, 1>;
57 using MatrixP = Eigen::Matrix<Real, N+M, N+M>;
58 using NewtonX = Optimist::RootFinder::Newton<Real, N>;
59 using NewtonK = Optimist::RootFinder::Newton<Real, N*S>;
74 using Time = Eigen::Vector<Real, Eigen::Dynamic>;
79 mutable Eigen::FullPivLU<MatrixP>
m_lu;
139 bool is_erk()
const {
return this->m_tableau.
type == Type::ERK;}
145 bool is_irk()
const {
return this->m_tableau.
type == Type::IRK;}
175 std::string
name()
const {
return this->m_tableau.
name;}
242 this->m_system = std::make_shared<ImplicitWrapper<Real, N, M>>(F, JF_x, JF_x_dot, h, Jh_x, in_domain);
264 this->m_system = std::make_shared<ImplicitWrapper<Real, N, M>>(
name, F, JF_x, JF_x_dot, h, Jh_x, in_domain);
282 this->m_system = std::make_shared<ExplicitWrapper<Real, N, M>>(f, Jf_x, h, Jh_x, in_domain);
302 this->m_system = std::make_shared<ExplicitWrapper<Real, N, M>>(
name, f, Jf_x, h, Jh_x, in_domain);
322 this->m_system = std::make_shared<LinearWrapper<Real, N, M>>(E,
A,
b, h, Jh_x, in_domain);
344 this->m_system = std::make_shared<LinearWrapper<Real, N, M>>(
name, E,
A,
b, h, Jh_x, in_domain);
366 this->m_system = std::make_shared<SemiExplicitWrapper<Real, N, M>>(
A, TA_x,
b, Jb_x, h, Jh_x, in_domain);
390 this->m_system = std::make_shared<SemiExplicitWrapper<Real, N, M>>(
name,
A, TA_x,
b, Jb_x, h, Jh_x, in_domain);
409 void absolute_tolerance(Real
const t_absolute_tolerance) {this->m_absolute_tolerance = t_absolute_tolerance;}
421 void relative_tolerance(Real
const t_relative_tolerance) {this->m_relative_tolerance = t_relative_tolerance;}
433 void safety_factor(Real
const t_safety_factor) {this->m_safety_factor = t_safety_factor;}
445 void min_safety_factor(Real
const t_min_safety_factor) {this->m_min_safety_factor = t_min_safety_factor;}
457 void max_safety_factor(Real
const t_max_safety_factor) {this->m_max_safety_factor = t_max_safety_factor;}
469 void min_step(Real
const t_min_step) {this->m_min_step = t_min_step;}
493 void adaptive(
bool t_adaptive) {this->m_adaptive = t_adaptive;}
516 this->m_verbose = t_verbose;
517 this->m_newtonX.verbose_mode(t_verbose);
518 this->m_newtonK.verbose_mode(t_verbose);
541 void reverse(
bool t_reverse) {this->m_reverse = t_reverse;}
576 {this->m_projection_tolerance = t_projection_tolerance;}
589 {this->m_max_projection_iterations = t_max_projection_iterations;}
601 void projection(
bool t_projection) {this->m_projection = t_projection;}
649 Real desired_error{this->m_absolute_tolerance + this->m_relative_tolerance *
650 std::max(x.array().abs().maxCoeff(), x_e.array().abs().maxCoeff())};
651 Real truncation_error{(x - x_e).array().abs().maxCoeff()};
652 return h_k * std::min(this->m_max_safety_factor, std::max(this->m_min_safety_factor,
653 this->m_safety_factor * std::pow(desired_error/truncation_error,
654 1.0/std::max(this->m_tableau.
order, this->m_tableau.order_e))));
663 std::ostringstream os;
665 <<
"Runge-Kutta method:\t" << this->
name() << std::endl
666 <<
"\t- order:\t" << this->
order() << std::endl
667 <<
"\t- stages:\t" << this->
stages() << std::endl
669 switch (this->
type()) {
670 case Type::ERK: os <<
"explicit";
break;
671 case Type::IRK: os <<
"implicit";
break;
672 case Type::DIRK: os <<
"diagonally implicit";
break;
676 <<
"\t- embedded:\t" << this->
is_embedded() << std::endl;
678 os <<
"\t- system:\t" << this->m_system->name() << std::endl;
680 os <<
"\t- system:\t" <<
"none" << std::endl;
722 for (
Integer i{0}; i < S; ++i) {
723 x_node = x_old + K(all, seqN(0, i)) * this->m_tableau.
A(i, seqN(0, i)).transpose();
724 if (!this->m_reverse) {
725 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));
730 if (!K.allFinite()) {
return false;}
733 x_new = x_old + K * this->m_tableau.
b;
736 if (this->m_adaptive && this->m_tableau.
is_embedded) {
737 VectorN x_emb = x_old + K * this->m_tableau.
b_e;
768 VectorN x_node(x + K(all, seqN(0, s)) * this->m_tableau.
A(s, seqN(0, s)).transpose());
769 if (!this->m_reverse) {
770 fun = this->m_system->F(x_node, K.col(s)/h, t + h * this->m_tableau.c(s));
772 fun = this->m_system->F_reverse(x_node, K.col(s)/h, t + h * this->m_tableau.c(s));
812 VectorN x_node(x + K(all, seqN(0, s)) * this->m_tableau.
A(s, seqN(0, s)).transpose());
813 if (!this->m_reverse) {
814 jac = this->m_system->JF_x_dot(x_node, K.col(s)/h, t + h * this->m_tableau.c(s)) / h;
816 jac = this->m_system->JF_x_dot_reverse(x_node, K.col(s)/h, t + h * this->m_tableau.c(s)) / h;
838 VectorN K_ini(VectorN::Zero());
841 for (
Integer s{0}; s < S; ++s) {
842 if (this->m_newtonX.solve(
843 [
this, s, &K, &x_old, t_old, h_old](
VectorN const & K_fun,
VectorN & fun)
844 {K.col(s) = K_fun; this->erk_implicit_function(s, x_old, t_old, h_old, K, fun);},
845 [
this, s, &K, &x_old, t_old, h_old](
VectorN const & K_jac,
MatrixN & jac)
846 {K.col(s) = K_jac; this->erk_implicit_jacobian(s, x_old, t_old, h_old, K, jac);},
855 x_new = x_old + K * this->m_tableau.
b;
858 if (this->m_adaptive && this->m_tableau.
is_embedded) {
859 VectorN x_emb(x_old + K * this->m_tableau.
b_e);
901 MatrixK K_mat{K.reshaped(N, S)};
903 for (
Integer i{0}; i < S; ++i) {
904 x_node = x + K_mat * this->m_tableau.
A.row(i).transpose();
905 if (!this->m_reverse) {
906 fun_mat.col(i) = this->m_system->F(x_node, K_mat.col(i)/h, t + h * this->m_tableau.c(i));
908 fun_mat.col(i) = this->m_system->F_reverse(x_node, K_mat.col(i)/h, t + h * this->m_tableau.c(i));
911 fun = fun_mat.reshaped(N*S, 1);
958 MatrixK K_mat{K.reshaped(N, S)};
962 auto idx = seqN(0, N), jdx = seqN(0, N);
963 for (
Integer i{0}; i < S; ++i) {
964 t_node = t + h * this->m_tableau.
c(i);
965 x_node = x + K_mat * this->m_tableau.
A.row(i).transpose();
968 x_dot_node = K_mat.col(i) / h;
969 if (!this->m_reverse) {
970 JF_x = this->m_system->JF_x(x_node, x_dot_node, t_node);
971 JF_x_dot = this->m_system->JF_x_dot(x_node, x_dot_node, t_node);
973 JF_x = this->m_system->JF_x_reverse(x_node, x_dot_node, t_node);
974 JF_x_dot = this->m_system->JF_x_dot_reverse(x_node, x_dot_node, t_node);
979 for (
Integer j{0}; j < S; ++j) {
982 jac(idx, jdx) = this->m_tableau.
A(i,j) * JF_x + JF_x_dot / h;
984 jac(idx, jdx) = this->m_tableau.
A(i,j) * JF_x;
1006 VectorK K_ini(VectorK::Zero());
1009 if (!this->m_newtonK.solve(
1010 [
this, &x_old, t_old, h_old](
VectorK const & K_fun,
VectorK & fun)
1011 {this->irk_function(x_old, t_old, h_old, K_fun, fun);},
1012 [
this, &x_old, t_old, h_old](
VectorK const & K_jac,
MatrixJ & jac)
1013 {this->irk_jacobian(x_old, t_old, h_old, K_jac, jac);},
1018 x_new = x_old + K.reshaped(N, S) * this->m_tableau.
b;
1021 if (this->m_adaptive && this->m_tableau.
is_embedded) {
1022 VectorN x_emb(x_old + K.reshaped(N, S) * this->m_tableau.b_e);
1062 VectorN x_node(x + K(all, seqN(0, n+1)) * this->m_tableau.
A(n, seqN(0, n+1)).transpose());
1063 if (!this->m_reverse) {
1064 fun = this->m_system->F(x_node, K.col(n)/h, t + h * this->m_tableau.c(n));
1066 fun = this->m_system->F_reverse(x_node, K.col(n)/h, t + h * this->m_tableau.c(n));
1108 Real t_node{t + h * this->m_tableau.
c(n)};
1109 VectorN x_node(x + K(all, seqN(0, n+1)) * this->m_tableau.
A(n, seqN(0, n+1)).transpose());
1110 VectorN x_dot_node(K.col(n)/h);
1111 if (!this->m_reverse) {
1112 jac = this->m_tableau.
A(n,n) * this->m_system->JF_x(x_node, x_dot_node, t_node) +
1113 this->m_system->JF_x_dot(x_node, x_dot_node, t_node) / h;
1115 jac = this->m_tableau.
A(n,n) * this->m_system->JF_x_reverse(x_node, x_dot_node, t_node) +
1116 this->m_system->JF_x_dot_reverse(x_node, x_dot_node, t_node) / h;
1137 VectorN K_ini(VectorN::Zero());
1140 for (
Integer n{0}; n < S; ++n) {
1141 if (this->m_newtonX.solve(
1142 [
this, n, &K, &x_old, t_old, h_old](
VectorN const & K_fun,
VectorN &fun)
1143 {K.col(n) = K_fun; this->dirk_function(n, x_old, t_old, h_old, K, fun);},
1144 [
this, n, &K, &x_old, t_old, h_old](
VectorN const & K_jac,
MatrixN &jac)
1145 {K.col(n) = K_jac; this->dirk_jacobian(n, x_old, t_old, h_old, K, jac);},
1154 x_new = x_old + K * this->m_tableau.
b;
1157 if (this->m_adaptive && this->m_tableau.
is_embedded) {
1158 VectorN x_emb(x_old + K * this->m_tableau.
b_e);
1175 bool step(
VectorN const & x_old, Real
const t_old, Real
const h_old,
VectorN & x_new, Real & h_new)
const
1177 #define CMD "Sandals::RungeKutta::step(...): "
1179 SANDALS_ASSERT(this->m_system->in_domain(x_old, t_old),
CMD "in " << this->m_tableau.name <<
1180 " solver, at t = " << t_old <<
", x = " << x_old.transpose() <<
", system out of domain.");
1182 if (this->
is_erk() && this->m_system->is_explicit()) {
1184 }
else if (this->
is_erk() && this->m_system->is_implicit()) {
1187 return this->
dirk_step(x_old, t_old, h_old, x_new, h_new);
1189 return this->
irk_step(x_old, t_old, h_old, x_new, h_new);
1210 #define CMD "Sandals::RungeKutta::advance(...): "
1214 h_old <<
", expected > 0.");
1217 if (!this->
step(x_old, t_old, h_old, x_new, h_new))
1220 Real t_tmp{t_old}, h_tmp{h_old / Real(2.0)};
1223 Integer max_k{this->m_max_substeps * this->m_max_substeps}, k{2};
1227 if (this->
step(x_tmp, t_tmp, h_tmp, x_new, h_new_tmp)) {
1233 if (k > 0 && k < max_k) {
1237 h_tmp = Real(2.0) * h_tmp;
1238 if (this->m_verbose) {
1240 ", integration succedded disable one substepping layer.");
1246 SANDALS_ASSERT(std::isfinite(x_tmp.maxCoeff()),
CMD "in " << this->m_tableau.name <<
1247 " solver, at t = " << t_tmp <<
", ||x||_inf = inf, computation interrupted.");
1254 t_tmp <<
", integration failed with h = " << h_tmp <<
", aborting.");
1259 "at t = " << t_tmp <<
", integration failed, adding substepping layer.");}
1274 if (this->m_projection) {
1276 if (this->
project(x_new, t_old + h_new, x_projected)) {
1277 x_new = x_projected;
1303 sol.
resize(t_mesh.size());
1306 sol.
t(0) = t_mesh(0);
1308 sol.
h.col(0) = this->m_system->h(ics, t_mesh(0));
1311 if (this->m_step_callback) {this->
m_step_callback(0, ics, t_mesh(0));}
1315 VectorN x_old_step(ics), x_new_step(ics);
1316 Real t_step{t_mesh(0)}, h_step{t_mesh(1) - t_mesh(0)}, h_tmp_step{h_step}, h_new_step;
1317 bool mesh_point_bool, saturation_bool;
1321 if (!this->
advance(x_old_step, t_step, h_step, x_new_step, h_new_step)) {
return false;}
1329 if (this->m_adaptive && this->m_tableau.
is_embedded && !mesh_point_bool && saturation_bool) {
1330 h_tmp_step = h_new_step;
1331 h_step = t_mesh(
step+1) - t_step;
1333 h_step = h_new_step;
1337 if (!this->m_adaptive || mesh_point_bool) {
1341 h_step = h_tmp_step;
1344 sol.
t(
step) = t_step;
1345 sol.
x.col(
step) = x_new_step;
1346 sol.
h.col(
step) = this->m_system->h(x_new_step, t_step);
1349 if (this->m_step_callback) {
1354 if (std::abs(t_step - t_mesh(last)) <
SQRT_EPSILON) {
break;}
1357 x_old_step = x_new_step;
1379 #define CMD "Sandals::RungeKutta::adaptive_solve(...): "
1384 return this->
solve(t_mesh, ics, sol);
1385 }
else if (!this->m_adaptive) {
1387 return this->
solve(t_mesh, ics, sol);
1391 Real t_step{t_mesh(0)}, h_step{t_mesh(1) - t_mesh(0)}, h_new_step, scale{100.0};
1392 Real h_min{std::max(this->m_min_step, h_step/scale)}, h_max{scale*h_step};
1394 Integer safety_length{
static_cast<Integer>(std::ceil(std::abs(t_mesh(last) - t_mesh(0))/(2.0*h_min)))};
1395 sol.
resize(safety_length);
1397 sol.
resize(t_mesh.size());
1401 sol.
t(0) = t_mesh(0);
1403 sol.
h.col(0) = this->m_system->h(ics, t_mesh(0));
1406 if (this->m_step_callback) {this->
m_step_callback(0, ics, t_mesh(0));}
1410 VectorN x_old_step(ics), x_new_step(ics);
1414 this->
advance(x_old_step, t_step, h_step, x_new_step, h_new_step);
1420 if (this->m_adaptive && this->m_tableau.
is_embedded) {
1421 h_step = std::max(std::min(h_new_step, h_max), h_min);
1430 sol.
t(
step) = t_step;
1431 sol.
x.col(
step) = x_new_step;
1432 sol.
h.col(
step) = this->m_system->h(x_new_step, t_step);
1435 if (this->m_step_callback) {
1440 if (std::abs(t_step - t_mesh(last)) <
SQRT_EPSILON) {
break;}
1441 else if (t_step + h_step > t_mesh(last)) {h_step = t_mesh(last) - t_step;}
1444 x_old_step = x_new_step;
1466 #define CMD "Sandals::RungeKutta::project(...): "
1477 A.template block<N, N>(0, 0) = MatrixN::Identity();
1478 for (
Integer k{0}; k < this->m_max_projection_iterations; ++k) {
1488 h = this->m_system->h(x_projected, t);
1489 Jh_x = this->m_system->Jh_x(x_projected, t);
1492 if (h.norm() < this->m_projection_tolerance) {
return true;}
1495 A.template block<N, M>(0, N) = Jh_x.transpose();
1496 A.template block<M, N>(N, 0) = Jh_x;
1497 b.template head<N>() = x - x_projected;
1498 b.template tail<M>() = -h;
1501 this->m_lu.compute(
A);
1503 x_step = this->m_lu.solve(
b);
1506 if (x_step.norm() < this->m_projection_tolerance * this->m_projection_tolerance) {
return false;}
1509 x_projected.noalias() += x_step(Eigen::seqN(0, N));
1532 std::vector<Integer>
const & projected_invariants,
VectorN & x_projected)
const
1534 #define CMD "Sandals::RungeKutta::project_ics(...): "
1548 A.block(0, 0, X, X) = MatrixX::Identity(X+H, X+H);
1549 Eigen::FullPivLU<MatrixX> lu;
1550 for (
Integer k{0}; k < this->m_max_projection_iterations; ++k) {
1560 h = this->m_system->h(x_projected, t);
1561 Jh_x = this->m_system->Jh_x(x_projected, t);
1564 h = h(projected_invariants);
1565 Jh_x = Jh_x(projected_invariants, projected_equations);
1568 if (h.norm() < this->m_projection_tolerance) {
return true;}
1571 A.block(0, X, X, H) = Jh_x.transpose();
1572 A.block(X, 0, H, X) = Jh_x;
1573 b.head(X) = x(projected_equations) - x_projected(projected_equations);
1579 x_step = this->m_lu.solve(
b);
1582 if (x_step.norm() < this->m_projection_tolerance * this->m_projection_tolerance) {
return false;}
1585 x_projected(projected_equations).noalias() += x_step;
1607 #define CMD "Sandals::RungeKutta::estimate_order(...): "
1611 for (
Integer i{0}; i < static_cast<Integer>(t_mesh.size()); ++i) {
1615 CMD "expected the same initial time.");
1617 CMD "expected the same final time.");
1620 for (
Integer j{1}; j < static_cast<Integer>(t_mesh[i].size()); ++j) {
1622 CMD "expected a fixed step.");
1629 VectorX h_vec(t_mesh.size()), e_vec(t_mesh.size());
1630 for (
Integer i{0}; i < static_cast<Integer>(t_mesh.size()); ++i) {
1632 "the" << i <<
"-th time mesh.");
1633 sol_ana = sol(sol_num.
t);
1635 CMD "expected the same number of states in analytical solution.");
1637 CMD "expected the same number of steps in analytical solution.");
1638 h_vec(i) = std::abs(sol_num.
t(1) - sol_num.
t(0));
1639 e_vec(i) = (sol_ana - sol_num.
x).array().abs().maxCoeff();
1645 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:160
static const FunctionH DefaultH
Definition Explicit.hh:258
std::function< MatrixJF(VectorF const &, Real const)> FunctionJF
Definition Explicit.hh:253
static const FunctionID DefaultID
Definition Explicit.hh:260
std::function< VectorH(VectorF const &, Real const)> FunctionH
Definition Explicit.hh:254
std::function< bool(VectorF const &, Real const)> FunctionID
Definition Explicit.hh:256
static const FunctionJH DefaultJH
Definition Explicit.hh:259
std::function< MatrixJH(VectorF const &, Real const)> FunctionJH
Definition Explicit.hh:255
std::function< VectorF(VectorF const &, Real const)> FunctionF
Definition Explicit.hh:252
Eigen::Matrix< Real, N, N > MatrixJF
Definition Implicit.hh:49
Eigen::Vector< Real, N > VectorF
Definition Implicit.hh:48
Eigen::Vector< Real, M > VectorH
Definition Implicit.hh:50
std::shared_ptr< Implicit< Real, N, M > > Pointer
Definition Implicit.hh:47
Eigen::Matrix< Real, M, N > MatrixJH
Definition Implicit.hh:51
std::function< bool(VectorF const &, Real const)> FunctionID
Definition Implicit.hh:283
std::function< VectorF(VectorF const &, VectorF const &, Real const)> FunctionF
Definition Implicit.hh:279
std::function< MatrixJF(VectorF const &, VectorF const &, Real const)> FunctionJF
Definition Implicit.hh:280
std::function< VectorH(VectorF const &, Real const)> FunctionH
Definition Implicit.hh:281
std::function< MatrixJH(VectorF const &, Real const)> FunctionJH
Definition Implicit.hh:282
static const FunctionID DefaultID
Definition Implicit.hh:287
static const FunctionH DefaultH
Definition Implicit.hh:285
static const FunctionJH DefaultJH
Definition Implicit.hh:286
std::function< VectorB(Real const)> FunctionB
Definition Linear.hh:217
static const FunctionH DefaultH
Definition Linear.hh:222
std::function< MatrixE(Real const)> FunctionE
Definition Linear.hh:215
static const FunctionJH DefaultJH
Definition Linear.hh:223
static const FunctionID DefaultID
Definition Linear.hh:224
std::function< MatrixJH(VectorF const &, Real const)> FunctionJH
Definition Linear.hh:219
std::function< bool(VectorF const &, Real const)> FunctionID
Definition Linear.hh:220
std::function< MatrixA(Real const)> FunctionA
Definition Linear.hh:216
std::function< VectorH(VectorF const &, Real const)> FunctionH
Definition Linear.hh:218
Real & min_step()
Definition RungeKutta.hh:463
bool projection()
Definition RungeKutta.hh:595
bool adaptive_mode()
Definition RungeKutta.hh:487
bool is_erk() const
Definition RungeKutta.hh:139
typename Tableau< Real, S >::Matrix MatrixS
Definition RungeKutta.hh:61
Real estimate_order(std::vector< VectorX > const &t_mesh, VectorN const &ics, std::function< MatrixX(VectorX)> &sol) const
Definition RungeKutta.hh:1603
void enable_reverse_mode()
Definition RungeKutta.hh:546
void projection(bool t_projection)
Definition RungeKutta.hh:601
void dirk_function(Integer n, VectorN const &x, Real const t, Real const h, MatrixK const &K, VectorN &fun) const
Definition RungeKutta.hh:1057
typename Tableau< Real, S >::Vector VectorS
Definition RungeKutta.hh:60
Real & min_safety_factor()
Definition RungeKutta.hh:439
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:357
bool reverse_mode()
Definition RungeKutta.hh:535
typename Implicit< Real, N, M >::MatrixJF MatrixN
Definition RungeKutta.hh:63
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:763
Real & safety_factor()
Definition RungeKutta.hh:427
Real m_absolute_tolerance
Definition RungeKutta.hh:83
Real m_min_safety_factor
Definition RungeKutta.hh:86
RungeKutta & operator=(RungeKutta const &)=delete
typename Implicit< Real, N, M >::MatrixJH MatrixM
Definition RungeKutta.hh:65
Eigen::Matrix< Real, N+M, 1 > VectorP
Definition RungeKutta.hh:56
Tableau< Real, S > & tableau()
Definition RungeKutta.hh:157
System system()
Definition RungeKutta.hh:217
VectorS b_embedded() const
Definition RungeKutta.hh:205
RungeKutta(const RungeKutta &)=delete
FunctionSC step_callback()
Definition RungeKutta.hh:557
bool advance(VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new) const
Definition RungeKutta.hh:1208
typename Implicit< Real, N, M >::Pointer System
Definition RungeKutta.hh:72
void enable_projection()
Definition RungeKutta.hh:606
RungeKutta(Tableau< Real, S > const &t_tableau, System t_system)
Definition RungeKutta.hh:124
Real m_relative_tolerance
Definition RungeKutta.hh:84
bool irk_step(VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new) const
Definition RungeKutta.hh:1003
Integer & max_substeps()
Definition RungeKutta.hh:475
typename Implicit< Real, N, M >::VectorF VectorN
Definition RungeKutta.hh:62
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:335
RungeKutta(Tableau< Real, S > const &t_tableau)
Definition RungeKutta.hh:114
Eigen::Matrix< Real, N+M, N+M > MatrixP
Definition RungeKutta.hh:57
void adaptive(bool t_adaptive)
Definition RungeKutta.hh:493
Integer order() const
Definition RungeKutta.hh:181
void projection_tolerance(Real const t_projection_tolerance)
Definition RungeKutta.hh:575
Eigen::Matrix< Real, N *S, N *S > MatrixJ
Definition RungeKutta.hh:55
void irk_jacobian(VectorN const &x, Real const t, Real const h, VectorK const &K, MatrixJ &jac) const
Definition RungeKutta.hh:950
VectorS b() const
Definition RungeKutta.hh:199
typename Tableau< Real, S >::Type Type
Definition RungeKutta.hh:73
NewtonX m_newtonX
Definition RungeKutta.hh:77
bool is_dirk() const
Definition RungeKutta.hh:151
VectorS c() const
Definition RungeKutta.hh:211
bool step(VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new) const
Definition RungeKutta.hh:1175
bool erk_implicit_step(VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new) const
Definition RungeKutta.hh:833
void reverse(bool t_reverse)
Definition RungeKutta.hh:541
bool m_reverse
Definition RungeKutta.hh:92
bool m_adaptive
Definition RungeKutta.hh:90
void max_safety_factor(Real const t_max_safety_factor)
Definition RungeKutta.hh:457
Eigen::Vector< Real, Eigen::Dynamic > VectorX
Definition RungeKutta.hh:51
std::string name() const
Definition RungeKutta.hh:175
Real m_min_step
Definition RungeKutta.hh:88
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:234
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:380
void enable_adaptive_mode()
Definition RungeKutta.hh:498
Real projection_tolerance()
Definition RungeKutta.hh:569
void verbose_mode(bool t_verbose)
Definition RungeKutta.hh:515
std::function< void(Integer const, VectorX const &, Real const)> FunctionSC
Definition RungeKutta.hh:66
void absolute_tolerance(Real const t_absolute_tolerance)
Definition RungeKutta.hh:409
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:807
Real relative_tolerance()
Definition RungeKutta.hh:415
bool verbose_mode()
Definition RungeKutta.hh:509
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:314
const Real SQRT_EPSILON
Definition RungeKutta.hh:70
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:275
Eigen::FullPivLU< MatrixP > m_lu
Definition RungeKutta.hh:79
System m_system
Definition RungeKutta.hh:82
void min_safety_factor(Real const t_min_safety_factor)
Definition RungeKutta.hh:445
Tableau< Real, S > m_tableau
Definition RungeKutta.hh:81
bool solve(VectorX const &t_mesh, VectorN const &ics, Solution< Real, N, M > &sol) const
Definition RungeKutta.hh:1298
bool adaptive_solve(VectorX const &t_mesh, VectorN const &ics, Solution< Real, N, M > &sol) const
Definition RungeKutta.hh:1374
void relative_tolerance(Real const t_relative_tolerance)
Definition RungeKutta.hh:421
void max_substeps(Integer const t_max_substeps)
Definition RungeKutta.hh:481
void min_step(Real const t_min_step)
Definition RungeKutta.hh:469
Eigen::Matrix< Real, N, S > MatrixK
Definition RungeKutta.hh:54
Integer stages() const
Definition RungeKutta.hh:169
bool dirk_step(VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new) const
Definition RungeKutta.hh:1133
void enable_verbose_mode()
Definition RungeKutta.hh:524
Integer & max_projection_iterations()
Definition RungeKutta.hh:582
bool is_irk() const
Definition RungeKutta.hh:145
Optimist::RootFinder::Newton< Real, N > NewtonX
Definition RungeKutta.hh:58
void step_callback(FunctionSC const &t_step_callback)
Definition RungeKutta.hh:563
Real m_projection_tolerance
Definition RungeKutta.hh:95
FunctionSC m_step_callback
Definition RungeKutta.hh:93
Real estimate_step(VectorN const &x, VectorN const &x_e, Real const h_k) const
Definition RungeKutta.hh:647
void dirk_jacobian(Integer n, VectorN const &x, Real const t, Real const h, MatrixK const &K, MatrixN &jac) const
Definition RungeKutta.hh:1103
void system(System t_system)
Definition RungeKutta.hh:223
void disable_adaptive_mode()
Definition RungeKutta.hh:503
Tableau< Real, S > const & tableau() const
Definition RungeKutta.hh:163
Real m_max_safety_factor
Definition RungeKutta.hh:87
bool m_verbose
Definition RungeKutta.hh:91
Real m_safety_factor
Definition RungeKutta.hh:85
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:1531
void disable_reverse_mode()
Definition RungeKutta.hh:551
bool erk_explicit_step(VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new) const
Definition RungeKutta.hh:713
Eigen::Vector< Real, N *S > VectorK
Definition RungeKutta.hh:53
void info(std::ostream &os)
Definition RungeKutta.hh:689
Type type() const
Definition RungeKutta.hh:133
void disable_projection()
Definition RungeKutta.hh:611
Eigen::Matrix< Real, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Definition RungeKutta.hh:52
void max_projection_iterations(Integer const t_max_projection_iterations)
Definition RungeKutta.hh:588
void safety_factor(Real const t_safety_factor)
Definition RungeKutta.hh:433
Eigen::Vector< Real, Eigen::Dynamic > Time
Definition RungeKutta.hh:74
void disable_verbose_mode()
Definition RungeKutta.hh:529
Integer m_max_substeps
Definition RungeKutta.hh:89
bool is_embedded() const
Definition RungeKutta.hh:187
std::string info() const
Definition RungeKutta.hh:661
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:294
bool has_system()
Definition RungeKutta.hh:397
void irk_function(VectorN const &x, Real const t, Real const h, VectorK const &K, VectorK &fun) const
Definition RungeKutta.hh:898
typename Implicit< Real, N, M >::VectorH VectorM
Definition RungeKutta.hh:64
MatrixS A() const
Definition RungeKutta.hh:193
Real absolute_tolerance()
Definition RungeKutta.hh:403
Optimist::RootFinder::Newton< Real, N *S > NewtonK
Definition RungeKutta.hh:59
Integer m_max_projection_iterations
Definition RungeKutta.hh:96
bool project(VectorN const &x, Real const t, VectorN &x_projected) const
Definition RungeKutta.hh:1464
NewtonK m_newtonK
Definition RungeKutta.hh:78
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:255
Real & max_safety_factor()
Definition RungeKutta.hh:451
bool m_projection
Definition RungeKutta.hh:97
std::function< MatrixJB(VectorF const &, Real const)> FunctionJB
Definition SemiExplicit.hh:277
static const FunctionH DefaultH
Definition SemiExplicit.hh:282
std::function< bool(VectorF const &, Real const)> FunctionID
Definition SemiExplicit.hh:280
std::function< VectorH(VectorF const &, Real const)> FunctionH
Definition SemiExplicit.hh:278
std::function< VectorB(VectorF const &, Real const)> FunctionB
Definition SemiExplicit.hh:276
std::function< TensorTA(VectorF const &, Real const)> FunctionTA
Definition SemiExplicit.hh:275
std::function< MatrixA(VectorF const &, Real const)> FunctionA
Definition SemiExplicit.hh:274
static const FunctionJH DefaultJH
Definition SemiExplicit.hh:283
static const FunctionID DefaultID
Definition SemiExplicit.hh:284
std::function< MatrixJH(VectorF const &, Real const)> FunctionJH
Definition SemiExplicit.hh:279
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