Sandals  v0.0.0
A C++ library for ODEs/DAEs integration
Loading...
Searching...
No Matches
RungeKutta.hh
Go to the documentation of this file.
1/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
2 * Copyright (c) 2025, Davide Stocco and Enrico Bertolazzi. *
3 * *
4 * The Sandals project is distributed under the BSD 2-Clause License. *
5 * *
6 * Davide Stocco Enrico Bertolazzi *
7 * University of Trento University of Trento *
8 * e-mail: davide.stocco@unitn.it e-mail: enrico.bertolazzi@unitn.it *
9\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
10
11#pragma once
12
13#ifndef SANDALS_RUNGEKUTTA_HH
14#define SANDALS_RUNGEKUTTA_HH
15
16#include <optional>
17
18#include <Sandals.hh>
19
24
25#include <Sandals/Tableau.hh>
26#include <Sandals/Solution.hh>
27
28namespace Sandals {
29
30 /*\
31 | ____ _ __ _ _
32 | | _ \ _ _ _ __ __ _ ___| |/ / _| |_| |_ __ _
33 | | |_) | | | | '_ \ / _` |/ _ \ ' / | | | __| __/ _` |
34 | | _ <| |_| | | | | (_| | __/ . \ |_| | |_| || (_| |
35 | |_| \_\\__,_|_| |_|\__, |\___|_|\_\__,_|\__|\__\__,_|
36 | |___/
37 \*/
38
50 template <typename Real, Integer S, Integer N, Integer M>
52 {
53 using VectorX = Eigen::Vector<Real, Eigen::Dynamic>;
54 using MatrixX = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
55 using VectorK = Eigen::Vector<Real, N*S>;
56 using MatrixK = Eigen::Matrix<Real, N, S>;
57 using MatrixJK = Eigen::Matrix<Real, N*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>;
69 using FunctionSC = std::function<void(Integer const, VectorX const &, Real const)>;
70
71 public:
73 const Real SQRT_EPSILON{std::sqrt(EPSILON)}; \
74
76 using Type = typename Tableau<Real, S>::Type;
77 using Time = Eigen::Vector<Real, Eigen::Dynamic>;
78
79 private:
82 mutable Eigen::FullPivLU<MatrixP> m_lu;
83
88 Real m_safety_factor{0.9};
91 Real m_min_step{EPSILON_HIGH};
93 bool m_adaptive{true};
94 bool m_verbose{false};
95 bool m_reverse{false};
97
98 Real m_projection_tolerance{EPSILON_HIGH};
100 bool m_projection{true};
101
102 public:
106 RungeKutta(const RungeKutta &) = delete;
107
111 RungeKutta & operator=(RungeKutta const & ) = delete;
112
117 RungeKutta(Tableau<Real, S> const & t_tableau)
118 : m_tableau(t_tableau) {
119 this->verbose_mode(this->m_verbose);
120 }
121
127 RungeKutta(Tableau<Real, S> const & t_tableau, System t_system)
128 : m_tableau(t_tableau), m_system(t_system) {
129 this->verbose_mode(this->m_verbose);
130 }
131
136 Type type() const {return this->m_tableau.type;}
137
142 bool is_erk() const {return this->m_tableau.type == Type::ERK;}
143
148 bool is_irk() const {return this->m_tableau.type == Type::IRK;}
149
154 bool is_dirk() const {return this->m_tableau.type == Type::DIRK;}
155
160 Tableau<Real, S> & tableau() {return this->m_tableau;}
161
166 Tableau<Real, S> const & tableau() const {return this->m_tableau;}
167
172 Integer stages() const {return S;}
173
178 std::string name() const {return this->m_tableau.name;}
179
184 Integer order() const {return this->m_tableau.order;}
185
190 bool is_embedded() const {return this->m_tableau.is_embedded;}
191
196 MatrixS A() const {return this->m_tableau.A;}
197
202 VectorS b() const {return this->m_tableau.b;}
203
208 VectorS b_embedded() const {return this->m_tableau.b_e;}
209
214 VectorS c() const {return this->m_tableau.c;}
215
220 System system() {return this->m_system;}
221
226 void system(System t_system) {this->m_system = t_system;}
227
244 ) {
245 this->m_system = std::make_shared<ImplicitWrapper<Real, N, M>>(F, JF_x, JF_x_dot, h, Jh_x, in_domain);
246 }
247
259 std::string name,
266 ) {
267 this->m_system = std::make_shared<ImplicitWrapper<Real, N, M>>(name, F, JF_x, JF_x_dot, h, Jh_x, in_domain);
268 }
269
284 ) {
285 this->m_system = std::make_shared<ExplicitWrapper<Real, N, M>>(f, Jf_x, h, Jh_x, in_domain);
286 }
287
298 std::string name,
304 ) {
305 this->m_system = std::make_shared<ExplicitWrapper<Real, N, M>>(name, f, Jf_x, h, Jh_x, in_domain);
306 }
307
324 ) {
325 this->m_system = std::make_shared<LinearWrapper<Real, N, M>>(E, A, b, h, Jh_x, in_domain);
326 }
327
339 std::string name,
346 ) {
347 this->m_system = std::make_shared<LinearWrapper<Real, N, M>>(name, E, A, b, h, Jh_x, in_domain);
348 }
349
371
395
400 bool has_system() {return this->m_system != nullptr;}
401
406 Real absolute_tolerance() {return this->m_absolute_tolerance;}
407
412 void absolute_tolerance(Real const t_absolute_tolerance) {this->m_absolute_tolerance = t_absolute_tolerance;}
413
418 Real relative_tolerance() {return this->m_relative_tolerance;}
419
424 void relative_tolerance(Real const t_relative_tolerance) {this->m_relative_tolerance = t_relative_tolerance;}
425
430 Real & safety_factor() {return this->m_safety_factor;}
431
436 void safety_factor(Real const t_safety_factor) {this->m_safety_factor = t_safety_factor;}
437
442 Real & min_safety_factor() {return this->m_min_safety_factor;}
443
448 void min_safety_factor(Real const t_min_safety_factor) {this->m_min_safety_factor = t_min_safety_factor;}
449
454 Real & max_safety_factor() {return this->m_max_safety_factor;}
455
460 void max_safety_factor(Real const t_max_safety_factor) {this->m_max_safety_factor = t_max_safety_factor;}
461
466 Real & min_step() {return this->m_min_step;}
467
472 void min_step(Real const t_min_step) {this->m_min_step = t_min_step;}
473
478 Integer & max_substeps() {return this->m_max_substeps;}
479
484 void max_substeps(Integer const t_max_substeps) {this->m_max_substeps = t_max_substeps;}
485
490 bool adaptive_mode() {return this->m_adaptive;}
491
496 void adaptive(bool t_adaptive) {this->m_adaptive = t_adaptive;}
497
501 void enable_adaptive_mode() {this->m_adaptive = true;}
502
506 void disable_adaptive_mode() {this->m_adaptive = false;}
507
512 bool verbose_mode() {return this->m_verbose;}
513
518 void verbose_mode(bool t_verbose) {
519 this->m_verbose = t_verbose;
520 this->m_newtonX.verbose_mode(t_verbose);
521 this->m_newtonK.verbose_mode(t_verbose);
522 }
523
527 void enable_verbose_mode() {this->verbose_mode(true);}
528
532 void disable_verbose_mode() {this->verbose_mode(false);}
533
538 bool reverse_mode() {return this->m_reverse;}
539
544 void reverse(bool t_reverse) {this->m_reverse = t_reverse;}
545
549 void enable_reverse_mode() {this->m_reverse = true;}
550
554 void disable_reverse_mode() {this->m_reverse = false;}
555
560 FunctionSC step_callback() {return this->m_step_callback;}
561
566 void step_callback(FunctionSC const & t_step_callback) {this->m_step_callback = t_step_callback;}
567
572 Real projection_tolerance() {return this->m_projection_tolerance;}
573
578 void projection_tolerance(Real const t_projection_tolerance)
579 {this->m_projection_tolerance = t_projection_tolerance;}
580
585 Integer & max_projection_iterations() {return this->m_max_projection_iterations;}
586
591 void max_projection_iterations(Integer const t_max_projection_iterations)
592 {this->m_max_projection_iterations = t_max_projection_iterations;}
593
598 bool projection() {return this->m_projection;}
599
604 void projection(bool t_projection) {this->m_projection = t_projection;}
605
609 void enable_projection() {this->m_projection = true;}
610
614 void disable_projection() {this->m_projection = false;}
615
650 Real estimate_step(VectorN const & x, VectorN const & x_e, Real const h_k) const
651 {
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))));
658 }
659
664 std::string info() const
665 {
666 std::ostringstream os;
667 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
671 << "\t- type:\t";
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;
676 }
677 os
678 << std::endl
679 << "\t- embedded:\t" << this->is_embedded() << std::endl;
680 if (this->has_system()) {
681 os << "\t- system:\t" << this->m_system->name() << std::endl;
682 } else {
683 os << "\t- system:\t" << "none" << std::endl;
684 }
685 return os.str();
686 }
687
692 void info(std::ostream &os) {os << this->info();}
693
694 /*\
695 | _____ ____ _ __
696 | | ____| _ \| |/ /
697 | | _| | |_) | ' /
698 | | |___| _ <| . \
699 | |_____|_| \_\_|\_\
700 |
701 \*/
702
717 bool erk_explicit_step(VectorN const & x_old, Real const t_old, Real const h_old, VectorN & x_new,
718 Real & h_new, MatrixK & K) const
719 {
720 using Eigen::all;
721 using Eigen::seqN;
722
723 // Compute the K variables in the case of an explicit method and explicit system
724 VectorN x_node;
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));
729 } else {
730 K.col(i) = h_old * static_cast<Explicit<Real, N, M> const *>(this->m_system.get())->f_reverse(x_node, t_old + h_old*this->m_tableau.c(i));
731 }
732 }
733 if (!K.allFinite()) {return false;}
734
735 // Perform the step and obtain the next state
736 x_new = x_old + K * this->m_tableau.b;
737
738 // Adapt next step
739 if (this->m_adaptive && this->m_tableau.is_embedded) {
740 VectorN x_emb = x_old + K * this->m_tableau.b_e;
741 h_new = this->estimate_step(x_new, x_emb, h_old);
742 }
743 return true;
744 }
745
766 void erk_implicit_function(Integer const s, VectorN const & x, Real const t, Real const h,
767 MatrixK const & K, VectorN & fun) const
768 {
769 using Eigen::all;
770 using Eigen::seqN;
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));
774 } else {
775 fun = this->m_system->F_reverse(x_node, K.col(s)/h, t + h * this->m_tableau.c(s));
776 }
777 }
778
810 void erk_implicit_jacobian(Integer const s, VectorN const & x, Real const t, Real const h,
811 MatrixK const & K, MatrixN & jac) const
812 {
813 using Eigen::all;
814 using Eigen::seqN;
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;
818 } else {
819 jac = this->m_system->JF_x_dot_reverse(x_node, K.col(s)/h, t + h * this->m_tableau.c(s)) / h;
820 }
821 }
822
837 bool erk_implicit_step(VectorN const & x_old, Real const t_old, Real const h_old, VectorN & x_new,
838 Real & h_new, MatrixK & K) const
839 {
840 VectorN K_sol;
841 VectorN K_ini(VectorN::Zero());
842
843 // Check if the solver converged
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);},
850 K_ini, K_sol)) {
851 K.col(s) = K_sol;
852 } else {
853 return false;
854 }
855 }
856
857 // Perform the step and obtain the next state
858 x_new = x_old + K * this->m_tableau.b;
859
860 // Adapt next step
861 if (this->m_adaptive && this->m_tableau.is_embedded) {
862 VectorN x_emb(x_old + K * this->m_tableau.b_e);
863 h_new = this->estimate_step(x_new, x_emb, h_old);
864 }
865 return true;
866 }
867
890 // void irk_jacobian_K(MatrixK const & K, MatrixJK & jac) const
891 // {
892//
893 // }
894
895 /*\
896 | ___ ____ _ __
897 | |_ _| _ \| |/ /
898 | | || |_) | ' /
899 | | || _ <| . \
900 | |___|_| \_\_|\_\
901 |
902 \*/
903
928 void irk_function(VectorN const & x, Real const t, Real const h, VectorK const & K, VectorK & fun) const
929 {
930 VectorN x_node;
931 MatrixK K_mat{K.reshaped(N, S)};
932 MatrixK fun_mat;
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));
937 } else {
938 fun_mat.col(i) = this->m_system->F_reverse(x_node, K_mat.col(i)/h, t + h * this->m_tableau.c(i));
939 }
940 }
941 fun = fun_mat.reshaped(N*S, 1);
942 }
943
980 void irk_jacobian(VectorN const & x, Real const t, Real const h, VectorK const & K, MatrixJ & jac) const
981 {
982 using Eigen::seqN;
983
984 // Reset the Jacobian matrix
985 jac.setZero();
986
987 // Loop through each equation of the system
988 MatrixK K_mat{K.reshaped(N, S)};
989 Real t_node;
990 VectorN x_node, x_dot_node;
991 MatrixN JF_x, JF_x_dot;
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();
996
997 // Compute the Jacobians with respect to x and x_dot
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);
1002 } else {
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);
1005 }
1006
1007 // Combine the Jacobians with respect to x and x_dot to obtain the Jacobian with respect to K
1008 idx = seqN(i*N, N);
1009 for (Integer j{0}; j < S; ++j) {
1010 jdx = seqN(j*N, N);
1011 if (i == j) {
1012 jac(idx, jdx) = this->m_tableau.A(i,j) * JF_x + JF_x_dot / h;
1013 } else {
1014 jac(idx, jdx) = this->m_tableau.A(i,j) * JF_x;
1015 }
1016 }
1017 }
1018 }
1019
1034 bool irk_step(VectorN const & x_old, Real const t_old, Real const h_old, VectorN & x_new,
1035 Real & h_new, MatrixK & K) const
1036 {
1037 VectorK K_vec;
1038 VectorK K_ini(VectorK::Zero());
1039
1040 // Check if the solver converged
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);},
1046 K_ini, K_vec))
1047 {return false;}
1048
1049 // Reshape the K vector to a matrix
1050 K = K_vec.reshaped(N, S);
1051 if (!K.allFinite()) {return false;}
1052
1053 // Perform the step and obtain the next state
1054 x_new = x_old + K * this->m_tableau.b;
1055
1056 // Adapt next step
1057 if (this->m_adaptive && this->m_tableau.is_embedded) {
1058 VectorN x_emb(x_old + K * this->m_tableau.b_e);
1059 h_new = this->estimate_step(x_new, x_emb, h_old);
1060 }
1061 return true;
1062 }
1063
1064 /*\
1065 | ____ ___ ____ _ __
1066 | | _ \_ _| _ \| |/ /
1067 | | | | | || |_) | ' /
1068 | | |_| | || _ <| . \
1069 | |____/___|_| \_\_|\_\
1070 |
1071 \*/
1072
1093 void dirk_function(Integer n, VectorN const & x, Real const t, Real const h, MatrixK const & K,
1094 VectorN & fun) const
1095 {
1096 using Eigen::all;
1097 using Eigen::seqN;
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));
1101 } else {
1102 fun = this->m_system->F_reverse(x_node, K.col(n)/h, t + h * this->m_tableau.c(n));
1103 }
1104 }
1105
1139 void dirk_jacobian(Integer n, VectorN const & x, Real const t, Real const h, MatrixK const & K,
1140 MatrixN & jac) const
1141 {
1142 using Eigen::all;
1143 using Eigen::seqN;
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;
1150 } else {
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;
1153 }
1154 }
1155
1170 bool dirk_step(VectorN const & x_old, Real const t_old, Real const h_old, VectorN & x_new,
1171 Real & h_new, MatrixK & K) const
1172 {
1173 VectorN K_sol;
1174 VectorN K_ini(VectorN::Zero());
1175
1176 // Check if the solver converged at each step
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);},
1183 K_ini, K_sol)) {
1184 K.col(n) = K_sol;
1185 } else {
1186 return false;
1187 }
1188 }
1189
1190 // Perform the step and obtain the next state
1191 x_new = x_old + K * this->m_tableau.b;
1192
1193 // Adapt next step
1194 if (this->m_adaptive && this->m_tableau.is_embedded) {
1195 VectorN x_emb(x_old + K * this->m_tableau.b_e);
1196 h_new = this->estimate_step(x_new, x_emb, h_old);
1197 }
1198 return true;
1199 }
1200
1213 bool step(VectorN const & x_old, Real const t_old, Real const h_old, VectorN & x_new, Real & h_new,
1214 MatrixK & K) const
1215 {
1216 #define CMD "Sandals::RungeKutta::step(...): "
1217
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.");
1220
1221 if (this->is_erk() && this->m_system->is_explicit()) {
1222 return this->erk_explicit_step(x_old, t_old, h_old, x_new, h_new, K);
1223 } else if (this->is_erk() && this->m_system->is_implicit()) {
1224 return this->erk_implicit_step(x_old, t_old, h_old, x_new, h_new, K);
1225 } else if (this->is_dirk()) {
1226 return this->dirk_step(x_old, t_old, h_old, x_new, h_new, K);
1227 } else {
1228 return this->irk_step(x_old, t_old, h_old, x_new, h_new, K);
1229 }
1230
1231 #undef CMD
1232 }
1233
1247 bool advance(VectorN const & x_old, Real const t_old, Real const h_old, VectorN & x_new, Real & h_new) const
1248 {
1249 #define CMD "Sandals::RungeKutta::advance(...): "
1250
1251 // Check step size
1252 SANDALS_ASSERT(h_old > Real(0.0), CMD "in " << this->m_tableau.name << " solver, h = "<<
1253 h_old << ", expected > 0.");
1254
1255 // If the integration step failed, try again with substepping
1256 MatrixK K;
1257 if (!this->step(x_old, t_old, h_old, x_new, h_new, K))
1258 {
1259 VectorN x_tmp(x_old);
1260 Real t_tmp{t_old}, h_tmp{h_old / Real(2.0)};
1261
1262 // Substepping logic
1263 Integer max_k{this->m_max_substeps * this->m_max_substeps}, k{2};
1264 Real h_new_tmp;
1265 while (k > 0) {
1266 // Calculate the next step with substepping logic
1267 if (this->step(x_tmp, t_tmp, h_tmp, x_new, h_new_tmp, K)) {
1268
1269 // Accept the step
1270 h_tmp = h_new_tmp;
1271
1272 // If substepping is enabled, double the step size
1273 if (k > 0 && k < max_k) {
1274 k -= 1;
1275 // If the substepping index is even, double the step size
1276 if (k % 2 == 0) {
1277 h_tmp = Real(2.0) * h_tmp;
1278 if (this->m_verbose) {
1279 SANDALS_WARNING(CMD "in " << this->m_tableau.name << " solver, at t = " << t_tmp <<
1280 ", integration succedded disable one substepping layer.");
1281 }
1282 }
1283 }
1284
1285 // Check the infinity norm of the solution
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.");
1288
1289 } else {
1290
1291 // If the substepping index is too high, abort the integration
1292 k += 2;
1293 SANDALS_ASSERT(k < max_k, CMD "in " << this->m_tableau.name << " solver, at t = " <<
1294 t_tmp << ", integration failed with h = " << h_tmp << ", aborting.");
1295 return false;
1296
1297 // Otherwise, try again with a smaller step
1298 if (this->m_verbose) {SANDALS_WARNING(CMD "in " << this->m_tableau.name << " solver, " <<
1299 "at t = " << t_tmp << ", integration failed, adding substepping layer.");}
1300 h_tmp /= Real(2.0);
1301 continue;
1302 }
1303
1304 // Store independent variable (or time)
1305 t_tmp += h_tmp;
1306 }
1307
1308 // Store output states substepping solutions
1309 x_new = x_tmp;
1310 h_new = h_tmp;
1311 }
1312
1313 // Project intermediate solution on the invariants
1314 if (this->m_projection) {
1315 VectorN x_projected;
1316 if (this->project(x_new, t_old + h_new, x_projected)) {
1317 x_new = x_projected;
1318 } else {
1319 return false;
1320 }
1321 }
1322 return true;
1323
1324 #undef CMD
1325 }
1326
1338 bool solve(VectorX const & t_mesh, VectorN const & ics, Solution<Real, N, M> & sol) const
1339 {
1340 using Eigen::last;
1341
1342 // Instantiate output
1343 sol.resize(t_mesh.size());
1344
1345 // Store initial conditions
1346 sol.t(0) = t_mesh(0);
1347 sol.x.col(0) = ics;
1348 sol.h.col(0) = this->m_system->h(ics, t_mesh(0));
1349
1350 // Callback on initial conditions
1351 if (this->m_step_callback) {this->m_step_callback(0, ics, t_mesh(0));}
1352
1353 // Update the current step
1354 Integer step{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;
1358
1359 while (true) {
1360 // Integrate system
1361 if (!this->advance(x_old_step, t_step, h_step, x_new_step, h_new_step)) {return false;}
1362
1363 // Update the current step
1364 t_step += h_step;
1365
1366 // Saturate the suggested timestep
1367 mesh_point_bool = std::abs(t_step - t_mesh(step+1)) < SQRT_EPSILON;
1368 saturation_bool = t_step + h_new_step > t_mesh(step+1) + SQRT_EPSILON;
1369 if (this->m_adaptive && this->m_tableau.is_embedded && !mesh_point_bool && saturation_bool) {
1370 h_tmp_step = h_new_step; // Used to store the previous step and keep the integration pace
1371 h_step = t_mesh(step+1) - t_step;
1372 } else {
1373 h_step = h_new_step;
1374 }
1375
1376 // Store solution if the step is a mesh point
1377 if (!this->m_adaptive || mesh_point_bool) {
1378
1379 // Update temporaries
1380 step += 1;
1381 h_step = h_tmp_step;
1382
1383 // Update outputs
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);
1387
1388 // Callback after the step is completed
1389 if (this->m_step_callback) {this->m_step_callback(step, x_new_step, t_step);}
1390
1391 // Check if the current step is the last one
1392 if (std::abs(t_step - t_mesh(last)) < SQRT_EPSILON) {break;}
1393
1394 // Update the previous step
1395 x_old_step = x_new_step;
1396 }
1397 }
1398 return true;
1399 }
1400
1412 bool adaptive_solve(VectorX const & t_mesh, VectorN const & ics, Solution<Real, N, M> & sol) const
1413 {
1414 using Eigen::all;
1415 using Eigen::last;
1416
1417 #define CMD "Sandals::RungeKutta::adaptive_solve(...): "
1418
1419 // Check if the adaptive method is enabled and the method is embedded
1420 if (!this->is_embedded()) {
1421 SANDALS_WARNING(CMD "the method is not embedded, using solve(...) method.");
1422 return this->solve(t_mesh, ics, sol);
1423 } else if (!this->m_adaptive) {
1424 SANDALS_WARNING(CMD "adaptive method is disabled, using solve(...) method.");
1425 return this->solve(t_mesh, ics, sol);
1426 }
1427
1428 // Instantiate output
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};
1431 if (this->m_tableau.is_embedded) {
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);
1434 } else {
1435 sol.resize(t_mesh.size());
1436 }
1437
1438 // Store initial conditions
1439 sol.t(0) = t_mesh(0);
1440 sol.x.col(0) = ics;
1441 sol.h.col(0) = this->m_system->h(ics, t_mesh(0));
1442
1443 // Callback on initial conditions
1444 if (this->m_step_callback) {this->m_step_callback(0, ics, t_mesh(0));}
1445
1446 // Instantiate temporary variables
1447 Integer step{0};
1448 VectorN x_old_step(ics), x_new_step(ics);
1449
1450 while (true) {
1451 // Integrate system
1452 this->advance(x_old_step, t_step, h_step, x_new_step, h_new_step);
1453
1454 // Update the current step
1455 t_step += h_step;
1456
1457 // Saturate the suggested timestep
1458 if (this->m_adaptive && this->m_tableau.is_embedded) {
1459 h_step = std::max(std::min(h_new_step, h_max), h_min);
1460 }
1461
1462 SANDALS_ASSERT(step < sol.size(), CMD "safety length exceeded.");
1463
1464 // Update temporaries
1465 step += 1;
1466
1467 // Update outputs
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);
1471
1472 // Callback after the step is completed
1473 if (this->m_step_callback) {this->m_step_callback(step, x_new_step, t_step);}
1474
1475 // Check if the current step is the last one
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;}
1478
1479 // Update the previous step
1480 x_old_step = x_new_step;
1481 }
1482
1483 // Resize the output
1485
1486 return true;
1487
1488 #undef CMD
1489 }
1490
1500 bool project(VectorN const & x, Real const t, VectorN & x_projected) const
1501 {
1502 #define CMD "Sandals::RungeKutta::project(...): "
1503
1504 // Check if there are any constraints
1505 x_projected = x;
1506 if (M > Integer(0)) {
1507
1508 VectorM h;
1509 MatrixM Jh_x;
1510 VectorP b, x_step;
1511 MatrixP A;
1512 A.setZero();
1513 A.template block<N, N>(0, 0) = MatrixN::Identity();
1514 for (Integer k{0}; k < this->m_max_projection_iterations; ++k) {
1515
1516 /* Standard projection method
1517 [A] {x} = {b}
1518 / I Jh_x^T \ / dx \ / x_t - x_k \
1519 | | | | = | |
1520 \ Jh_x 0 / \ lambda / \ -h /
1521 */
1522
1523 // Evaluate the invariants vector and its Jacobian
1524 h = this->m_system->h(x_projected, t);
1525 Jh_x = this->m_system->Jh_x(x_projected, t);
1526
1527 // Check if the solution is found
1528 if (h.norm() < this->m_projection_tolerance) {return true;}
1529
1530 // Build the linear system
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;
1535
1536 // Compute the solution of the linear system
1537 this->m_lu.compute(A);
1538 SANDALS_ASSERT(this->m_lu.rank() == N+M, CMD "singular Jacobian detected.");
1539 x_step = this->m_lu.solve(b);
1540
1541 // Check if the step is too small
1542 if (x_step.norm() < this->m_projection_tolerance * this->m_projection_tolerance) {return false;}
1543
1544 // Update the solution
1545 x_projected.noalias() += x_step(Eigen::seqN(0, N));
1546 }
1547 if (this->m_verbose) {SANDALS_WARNING(CMD "maximum number of iterations reached.");}
1548 return false;
1549 } else {
1550 return true;
1551 }
1552
1553 #undef CMD
1554 }
1555
1567 bool project_ics(VectorN const & x, Real const t, std::vector<Integer> const & projected_equations,
1568 std::vector<Integer> const & projected_invariants, VectorN & x_projected) const
1569 {
1570 #define CMD "Sandals::RungeKutta::project_ics(...): "
1571
1572 Integer X{static_cast<Integer>(projected_equations.size())};
1573 Integer H{static_cast<Integer>(projected_invariants.size())};
1574
1575 // Check if there are any constraints
1576 x_projected = x;
1577 if (H > Integer(0)) {
1578
1579 VectorM h;
1580 MatrixM Jh_x;
1581 VectorX b(X+H), x_step(X+H);
1582 MatrixX A(X+H, X+H);
1583 A.setZero();
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) {
1587
1588 /* Standard projection method
1589 [A] {x} = {b}
1590 / I Jh_x^T \ / dx \ / x_t - x_k \
1591 | | | | = | |
1592 \ Jh_x 0 / \ lambda / \ -h /
1593 */
1594
1595 // Evaluate the invariants vector and its Jacobian
1596 h = this->m_system->h(x_projected, t);
1597 Jh_x = this->m_system->Jh_x(x_projected, t);
1598
1599 // Select only the projected invariants
1600 h = h(projected_invariants);
1601 Jh_x = Jh_x(projected_invariants, projected_equations);
1602
1603 // Check if the solution is found
1604 if (h.norm() < this->m_projection_tolerance) {return true;}
1605
1606 // Build the linear system
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);
1610 b.tail(H) = -h;
1611
1612 // Compute the solution of the linear system
1613 lu.compute(A);
1614 SANDALS_ASSERT(lu.rank() == X+H, CMD "singular Jacobian detected.");
1615 x_step = this->m_lu.solve(b);
1616
1617 // Check if the step is too small
1618 if (x_step.norm() < this->m_projection_tolerance * this->m_projection_tolerance) {return false;}
1619
1620 // Update the solution
1621 x_projected(projected_equations).noalias() += x_step;
1622 }
1623 if (this->m_verbose) {SANDALS_WARNING(CMD "maximum number of iterations reached.");}
1624 return false;
1625 } else {
1626 return true;
1627 }
1628
1629 #undef CMD
1630 }
1631
1639 Real estimate_order(std::vector<VectorX> const & t_mesh, VectorN const & ics, std::function<MatrixX(VectorX)> & sol) const
1640 {
1641 using Eigen::last;
1642
1643 #define CMD "Sandals::RungeKutta::estimate_order(...): "
1644
1645 SANDALS_ASSERT(t_mesh.size() > Integer(1), CMD "expected at least two time meshes.");
1646
1647 for (Integer i{0}; i < static_cast<Integer>(t_mesh.size()); ++i) {
1648
1649 // Check if the time mesh has the same initial and final time
1650 SANDALS_ASSERT((t_mesh[0](0) - t_mesh[i](0)) < SQRT_EPSILON,
1651 CMD "expected the same initial time.");
1652 SANDALS_ASSERT((t_mesh[0](last) - t_mesh[i](last)) < SQRT_EPSILON,
1653 CMD "expected the same final time.");
1654
1655 // Check if the mesh step is fixed all over the time mesh
1656 for (Integer j{1}; j < static_cast<Integer>(t_mesh[i].size()); ++j) {
1657 SANDALS_ASSERT((t_mesh[i](j) - t_mesh[i](j-1)) - (t_mesh[i](1) - t_mesh[i](0)) < SQRT_EPSILON,
1658 CMD "expected a fixed step.");
1659 }
1660 }
1661
1662 // Solve the system for each time scale
1663 Solution<Real, N, M> sol_num;
1664 MatrixX sol_ana;
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) {
1667 SANDALS_ASSERT(this->solve(t_mesh[i], ics, sol_num), CMD "failed to solve the system for " <<
1668 "the" << i << "-th time mesh.");
1669 sol_ana = sol(sol_num.t);
1670 SANDALS_ASSERT(sol_ana.rows() == sol_num.x.rows(),
1671 CMD "expected the same number of states in analytical solution.");
1672 SANDALS_ASSERT(sol_ana.cols() == sol_num.x.cols(),
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();
1676 }
1677
1678 // Compute the order of the method thorugh least squares
1679 VectorX A(h_vec.array().log());
1680 VectorX b(e_vec.array().log());
1681 return ((A.transpose() * A).ldlt().solve(A.transpose() * b))(0);
1682
1683 #undef CMD
1684 }
1685
1686 }; // class RungeKutta
1687
1688} // namespace Sandals
1689
1690#endif // SANDALS_RUNGEKUTTA_HH
#define CMD
#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