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 <Sandals.hh>
17
22
23#include <Sandals/Tableau.hh>
24#include <Sandals/Solution.hh>
25
26namespace Sandals {
27
28 /*\
29 | ____ _ __ _ _
30 | | _ \ _ _ _ __ __ _ ___| |/ / _| |_| |_ __ _
31 | | |_) | | | | '_ \ / _` |/ _ \ ' / | | | __| __/ _` |
32 | | _ <| |_| | | | | (_| | __/ . \ |_| | |_| || (_| |
33 | |_| \_\\__,_|_| |_|\__, |\___|_|\_\__,_|\__|\__\__,_|
34 | |___/
35 \*/
36
48 template <typename Real, Integer S, Integer N, Integer M>
50 {
51 using VectorX = Eigen::Vector<Real, Eigen::Dynamic>;
52 using MatrixX = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
53 using VectorK = Eigen::Vector<Real, N*S>;
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>;
66 using FunctionSC = std::function<void(Integer const, VectorX const &, Real const)>;
67
68 public:
70 const Real SQRT_EPSILON{std::sqrt(EPSILON)}; \
71
73 using Type = typename Tableau<Real, S>::Type;
74 using Time = Eigen::Vector<Real, Eigen::Dynamic>;
75
76 private:
79 mutable Eigen::FullPivLU<MatrixP> m_lu;
80
85 Real m_safety_factor{0.9};
88 Real m_min_step{EPSILON_HIGH};
90 bool m_adaptive{true};
91 bool m_verbose{false};
92 bool m_reverse{false};
94
95 Real m_projection_tolerance{EPSILON_HIGH};
97 bool m_projection{true};
98
99 public:
103 RungeKutta(const RungeKutta &) = delete;
104
108 RungeKutta & operator=(RungeKutta const & ) = delete;
109
114 RungeKutta(Tableau<Real, S> const & t_tableau)
115 : m_tableau(t_tableau) {
116 this->verbose_mode(this->m_verbose);
117 }
118
124 RungeKutta(Tableau<Real, S> const & t_tableau, System t_system)
125 : m_tableau(t_tableau), m_system(t_system) {
126 this->verbose_mode(this->m_verbose);
127 }
128
133 Type type() const {return this->m_tableau.type;}
134
139 bool is_erk() const {return this->m_tableau.type == Type::ERK;}
140
145 bool is_irk() const {return this->m_tableau.type == Type::IRK;}
146
151 bool is_dirk() const {return this->m_tableau.type == Type::DIRK;}
152
157 Tableau<Real, S> & tableau() {return this->m_tableau;}
158
163 Tableau<Real, S> const & tableau() const {return this->m_tableau;}
164
169 Integer stages() const {return S;}
170
175 std::string name() const {return this->m_tableau.name;}
176
181 Integer order() const {return this->m_tableau.order;}
182
187 bool is_embedded() const {return this->m_tableau.is_embedded;}
188
193 MatrixS A() const {return this->m_tableau.A;}
194
199 VectorS b() const {return this->m_tableau.b;}
200
205 VectorS b_embedded() const {return this->m_tableau.b_e;}
206
211 VectorS c() const {return this->m_tableau.c;}
212
217 System system() {return this->m_system;}
218
223 void system(System t_system) {this->m_system = t_system;}
224
241 ) {
242 this->m_system = std::make_shared<ImplicitWrapper<Real, N, M>>(F, JF_x, JF_x_dot, h, Jh_x, in_domain);
243 }
244
256 std::string name,
263 ) {
264 this->m_system = std::make_shared<ImplicitWrapper<Real, N, M>>(name, F, JF_x, JF_x_dot, h, Jh_x, in_domain);
265 }
266
281 ) {
282 this->m_system = std::make_shared<ExplicitWrapper<Real, N, M>>(f, Jf_x, h, Jh_x, in_domain);
283 }
284
295 std::string name,
301 ) {
302 this->m_system = std::make_shared<ExplicitWrapper<Real, N, M>>(name, f, Jf_x, h, Jh_x, in_domain);
303 }
304
321 ) {
322 this->m_system = std::make_shared<LinearWrapper<Real, N, M>>(E, A, b, h, Jh_x, in_domain);
323 }
324
336 std::string name,
343 ) {
344 this->m_system = std::make_shared<LinearWrapper<Real, N, M>>(name, E, A, b, h, Jh_x, in_domain);
345 }
346
368
392
397 bool has_system() {return this->m_system != nullptr;}
398
403 Real absolute_tolerance() {return this->m_absolute_tolerance;}
404
409 void absolute_tolerance(Real const t_absolute_tolerance) {this->m_absolute_tolerance = t_absolute_tolerance;}
410
415 Real relative_tolerance() {return this->m_relative_tolerance;}
416
421 void relative_tolerance(Real const t_relative_tolerance) {this->m_relative_tolerance = t_relative_tolerance;}
422
427 Real & safety_factor() {return this->m_safety_factor;}
428
433 void safety_factor(Real const t_safety_factor) {this->m_safety_factor = t_safety_factor;}
434
439 Real & min_safety_factor() {return this->m_min_safety_factor;}
440
445 void min_safety_factor(Real const t_min_safety_factor) {this->m_min_safety_factor = t_min_safety_factor;}
446
451 Real & max_safety_factor() {return this->m_max_safety_factor;}
452
457 void max_safety_factor(Real const t_max_safety_factor) {this->m_max_safety_factor = t_max_safety_factor;}
458
463 Real & min_step() {return this->m_min_step;}
464
469 void min_step(Real const t_min_step) {this->m_min_step = t_min_step;}
470
475 Integer & max_substeps() {return this->m_max_substeps;}
476
481 void max_substeps(Integer const t_max_substeps) {this->m_max_substeps = t_max_substeps;}
482
487 bool adaptive_mode() {return this->m_adaptive;}
488
493 void adaptive(bool t_adaptive) {this->m_adaptive = t_adaptive;}
494
498 void enable_adaptive_mode() {this->m_adaptive = true;}
499
503 void disable_adaptive_mode() {this->m_adaptive = false;}
504
509 bool verbose_mode() {return this->m_verbose;}
510
515 void verbose_mode(bool t_verbose) {
516 this->m_verbose = t_verbose;
517 this->m_newtonX.verbose_mode(t_verbose);
518 this->m_newtonK.verbose_mode(t_verbose);
519 }
520
524 void enable_verbose_mode() {this->verbose_mode(true);}
525
529 void disable_verbose_mode() {this->verbose_mode(false);}
530
535 bool reverse_mode() {return this->m_reverse;}
536
541 void reverse(bool t_reverse) {this->m_reverse = t_reverse;}
542
546 void enable_reverse_mode() {this->m_reverse = true;}
547
551 void disable_reverse_mode() {this->m_reverse = false;}
552
557 FunctionSC step_callback() {return this->m_step_callback;}
558
563 void step_callback(FunctionSC const & t_step_callback) {this->m_step_callback = t_step_callback;}
564
569 Real projection_tolerance() {return this->m_projection_tolerance;}
570
575 void projection_tolerance(Real const t_projection_tolerance)
576 {this->m_projection_tolerance = t_projection_tolerance;}
577
582 Integer & max_projection_iterations() {return this->m_max_projection_iterations;}
583
588 void max_projection_iterations(Integer const t_max_projection_iterations)
589 {this->m_max_projection_iterations = t_max_projection_iterations;}
590
595 bool projection() {return this->m_projection;}
596
601 void projection(bool t_projection) {this->m_projection = t_projection;}
602
606 void enable_projection() {this->m_projection = true;}
607
611 void disable_projection() {this->m_projection = false;}
612
647 Real estimate_step(VectorN const & x, VectorN const & x_e, Real const h_k) const
648 {
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))));
655 }
656
661 std::string info() const
662 {
663 std::ostringstream os;
664 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
668 << "\t- type:\t";
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;
673 }
674 os
675 << std::endl
676 << "\t- embedded:\t" << this->is_embedded() << std::endl;
677 if (this->has_system()) {
678 os << "\t- system:\t" << this->m_system->name() << std::endl;
679 } else {
680 os << "\t- system:\t" << "none" << std::endl;
681 }
682 return os.str();
683 }
684
689 void info(std::ostream &os) {os << this->info();}
690
691 /*\
692 | _____ ____ _ __
693 | | ____| _ \| |/ /
694 | | _| | |_) | ' /
695 | | |___| _ <| . \
696 | |_____|_| \_\_|\_\
697 |
698 \*/
699
713 bool erk_explicit_step(VectorN const & x_old, Real const t_old, Real const h_old, VectorN & x_new,
714 Real & h_new) const
715 {
716 using Eigen::all;
717 using Eigen::seqN;
718
719 // Compute the K variables in the case of an explicit method and explicit system
720 VectorN x_node;
721 MatrixK K;
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));
726 } else {
727 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));
728 }
729 }
730 if (!K.allFinite()) {return false;}
731
732 // Perform the step and obtain the next state
733 x_new = x_old + K * this->m_tableau.b;
734
735 // Adapt next step
736 if (this->m_adaptive && this->m_tableau.is_embedded) {
737 VectorN x_emb = x_old + K * this->m_tableau.b_e;
738 h_new = this->estimate_step(x_new, x_emb, h_old);
739 }
740 return true;
741 }
742
763 void erk_implicit_function(Integer const s, VectorN const & x, Real const t, Real const h,
764 MatrixK const & K, VectorN & fun) const
765 {
766 using Eigen::all;
767 using Eigen::seqN;
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));
771 } else {
772 fun = this->m_system->F_reverse(x_node, K.col(s)/h, t + h * this->m_tableau.c(s));
773 }
774 }
775
807 void erk_implicit_jacobian(Integer const s, VectorN const & x, Real const t, Real const h,
808 MatrixK const & K, MatrixN & jac) const
809 {
810 using Eigen::all;
811 using Eigen::seqN;
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;
815 } else {
816 jac = this->m_system->JF_x_dot_reverse(x_node, K.col(s)/h, t + h * this->m_tableau.c(s)) / h;
817 }
818 }
819
833 bool erk_implicit_step(VectorN const & x_old, Real const t_old, Real const h_old, VectorN & x_new,
834 Real & h_new) const
835 {
836 MatrixK K;
837 VectorN K_sol;
838 VectorN K_ini(VectorN::Zero());
839
840 // Check if the solver converged
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);},
847 K_ini, K_sol)) {
848 K.col(s) = K_sol;
849 } else {
850 return false;
851 }
852 }
853
854 // Perform the step and obtain the next state
855 x_new = x_old + K * this->m_tableau.b;
856
857 // Adapt next step
858 if (this->m_adaptive && this->m_tableau.is_embedded) {
859 VectorN x_emb(x_old + K * this->m_tableau.b_e);
860 h_new = this->estimate_step(x_new, x_emb, h_old);
861 }
862 return true;
863 }
864
865 /*\
866 | ___ ____ _ __
867 | |_ _| _ \| |/ /
868 | | || |_) | ' /
869 | | || _ <| . \
870 | |___|_| \_\_|\_\
871 |
872 \*/
873
898 void irk_function(VectorN const & x, Real const t, Real const h, VectorK const & K, VectorK & fun) const
899 {
900 VectorN x_node;
901 MatrixK K_mat{K.reshaped(N, S)};
902 MatrixK fun_mat;
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));
907 } else {
908 fun_mat.col(i) = this->m_system->F_reverse(x_node, K_mat.col(i)/h, t + h * this->m_tableau.c(i));
909 }
910 }
911 fun = fun_mat.reshaped(N*S, 1);
912 }
913
950 void irk_jacobian(VectorN const & x, Real const t, Real const h, VectorK const & K, MatrixJ & jac) const
951 {
952 using Eigen::seqN;
953
954 // Reset the Jacobian matrix
955 jac.setZero();
956
957 // Loop through each equation of the system
958 MatrixK K_mat{K.reshaped(N, S)};
959 Real t_node;
960 VectorN x_node, x_dot_node;
961 MatrixN JF_x, JF_x_dot;
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();
966
967 // Compute the Jacobians with respect to x and x_dot
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);
972 } else {
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);
975 }
976
977 // Combine the Jacobians with respect to x and x_dot to obtain the Jacobian with respect to K
978 idx = seqN(i*N, N);
979 for (Integer j{0}; j < S; ++j) {
980 jdx = seqN(j*N, N);
981 if (i == j) {
982 jac(idx, jdx) = this->m_tableau.A(i,j) * JF_x + JF_x_dot / h;
983 } else {
984 jac(idx, jdx) = this->m_tableau.A(i,j) * JF_x;
985 }
986 }
987 }
988 }
989
1003 bool irk_step(VectorN const & x_old, Real const t_old, Real const h_old, VectorN & x_new, Real & h_new) const
1004 {
1005 VectorK K;
1006 VectorK K_ini(VectorK::Zero());
1007
1008 // Check if the solver converged
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);},
1014 K_ini, K))
1015 {return false;}
1016
1017 // Perform the step and obtain the next state
1018 x_new = x_old + K.reshaped(N, S) * this->m_tableau.b;
1019
1020 // Adapt next step
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);
1023 h_new = this->estimate_step(x_new, x_emb, h_old);
1024 }
1025 return true;
1026 }
1027
1028 /*\
1029 | ____ ___ ____ _ __
1030 | | _ \_ _| _ \| |/ /
1031 | | | | | || |_) | ' /
1032 | | |_| | || _ <| . \
1033 | |____/___|_| \_\_|\_\
1034 |
1035 \*/
1036
1057 void dirk_function(Integer n, VectorN const & x, Real const t, Real const h, MatrixK const & K,
1058 VectorN & fun) const
1059 {
1060 using Eigen::all;
1061 using Eigen::seqN;
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));
1065 } else {
1066 fun = this->m_system->F_reverse(x_node, K.col(n)/h, t + h * this->m_tableau.c(n));
1067 }
1068 }
1069
1103 void dirk_jacobian(Integer n, VectorN const & x, Real const t, Real const h, MatrixK const & K,
1104 MatrixN & jac) const
1105 {
1106 using Eigen::all;
1107 using Eigen::seqN;
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;
1114 } else {
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;
1117 }
1118 }
1119
1133 bool dirk_step(VectorN const & x_old, Real const t_old, Real const h_old, VectorN & x_new, Real & h_new) const
1134 {
1135 MatrixK K;
1136 VectorN K_sol;
1137 VectorN K_ini(VectorN::Zero());
1138
1139 // Check if the solver converged at each step
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);},
1146 K_ini, K_sol)) {
1147 K.col(n) = K_sol;
1148 } else {
1149 return false;
1150 }
1151 }
1152
1153 // Perform the step and obtain the next state
1154 x_new = x_old + K * this->m_tableau.b;
1155
1156 // Adapt next step
1157 if (this->m_adaptive && this->m_tableau.is_embedded) {
1158 VectorN x_emb(x_old + K * this->m_tableau.b_e);
1159 h_new = this->estimate_step(x_new, x_emb, h_old);
1160 }
1161 return true;
1162 }
1163
1175 bool step(VectorN const & x_old, Real const t_old, Real const h_old, VectorN & x_new, Real & h_new) const
1176 {
1177 #define CMD "Sandals::RungeKutta::step(...): "
1178
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.");
1181
1182 if (this->is_erk() && this->m_system->is_explicit()) {
1183 return this->erk_explicit_step(x_old, t_old, h_old, x_new, h_new);
1184 } else if (this->is_erk() && this->m_system->is_implicit()) {
1185 return this->erk_implicit_step(x_old, t_old, h_old, x_new, h_new);
1186 } else if (this->is_dirk()) {
1187 return this->dirk_step(x_old, t_old, h_old, x_new, h_new);
1188 } else {
1189 return this->irk_step(x_old, t_old, h_old, x_new, h_new);
1190 }
1191
1192 #undef CMD
1193 }
1194
1208 bool advance(VectorN const & x_old, Real const t_old, Real const h_old, VectorN & x_new, Real & h_new) const
1209 {
1210 #define CMD "Sandals::RungeKutta::advance(...): "
1211
1212 // Check step size
1213 SANDALS_ASSERT(h_old > Real(0.0), CMD "in " << this->m_tableau.name << " solver, h = "<<
1214 h_old << ", expected > 0.");
1215
1216 // If the integration step failed, try again with substepping
1217 if (!this->step(x_old, t_old, h_old, x_new, h_new))
1218 {
1219 VectorN x_tmp(x_old);
1220 Real t_tmp{t_old}, h_tmp{h_old / Real(2.0)};
1221
1222 // Substepping logic
1223 Integer max_k{this->m_max_substeps * this->m_max_substeps}, k{2};
1224 Real h_new_tmp;
1225 while (k > 0) {
1226 // Calculate the next step with substepping logic
1227 if (this->step(x_tmp, t_tmp, h_tmp, x_new, h_new_tmp)) {
1228
1229 // Accept the step
1230 h_tmp = h_new_tmp;
1231
1232 // If substepping is enabled, double the step size
1233 if (k > 0 && k < max_k) {
1234 k -= 1;
1235 // If the substepping index is even, double the step size
1236 if (k % 2 == 0) {
1237 h_tmp = Real(2.0) * h_tmp;
1238 if (this->m_verbose) {
1239 SANDALS_WARNING(CMD "in " << this->m_tableau.name << " solver, at t = " << t_tmp <<
1240 ", integration succedded disable one substepping layer.");
1241 }
1242 }
1243 }
1244
1245 // Check the infinity norm of the solution
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.");
1248
1249 } else {
1250
1251 // If the substepping index is too high, abort the integration
1252 k += 2;
1253 SANDALS_ASSERT(k < max_k, CMD "in " << this->m_tableau.name << " solver, at t = " <<
1254 t_tmp << ", integration failed with h = " << h_tmp << ", aborting.");
1255 return false;
1256
1257 // Otherwise, try again with a smaller step
1258 if (this->m_verbose) {SANDALS_WARNING(CMD "in " << this->m_tableau.name << " solver, " <<
1259 "at t = " << t_tmp << ", integration failed, adding substepping layer.");}
1260 h_tmp /= Real(2.0);
1261 continue;
1262 }
1263
1264 // Store independent variable (or time)
1265 t_tmp += h_tmp;
1266 }
1267
1268 // Store output states substepping solutions
1269 x_new = x_tmp;
1270 h_new = h_tmp;
1271 }
1272
1273 // Project intermediate solution on the invariants
1274 if (this->m_projection) {
1275 VectorN x_projected;
1276 if (this->project(x_new, t_old + h_new, x_projected)) {
1277 x_new = x_projected;
1278 } else {
1279 return false;
1280 }
1281 }
1282 return true;
1283
1284 #undef CMD
1285 }
1286
1298 bool solve(VectorX const & t_mesh, VectorN const & ics, Solution<Real, N, M> & sol) const
1299 {
1300 using Eigen::last;
1301
1302 // Instantiate output
1303 sol.resize(t_mesh.size());
1304
1305 // Store initial conditions
1306 sol.t(0) = t_mesh(0);
1307 sol.x.col(0) = ics;
1308 sol.h.col(0) = this->m_system->h(ics, t_mesh(0));
1309
1310 // Callback on initial conditions
1311 if (this->m_step_callback) {this->m_step_callback(0, ics, t_mesh(0));}
1312
1313 // Update the current step
1314 Integer step{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;
1318
1319 while (true) {
1320 // Integrate system
1321 if (!this->advance(x_old_step, t_step, h_step, x_new_step, h_new_step)) {return false;}
1322
1323 // Update the current step
1324 t_step += h_step;
1325
1326 // Saturate the suggested timestep
1327 mesh_point_bool = std::abs(t_step - t_mesh(step+1)) < SQRT_EPSILON;
1328 saturation_bool = t_step + h_new_step > t_mesh(step+1) + SQRT_EPSILON;
1329 if (this->m_adaptive && this->m_tableau.is_embedded && !mesh_point_bool && saturation_bool) {
1330 h_tmp_step = h_new_step; // Used to store the previous step and keep the integration pace
1331 h_step = t_mesh(step+1) - t_step;
1332 } else {
1333 h_step = h_new_step;
1334 }
1335
1336 // Store solution if the step is a mesh point
1337 if (!this->m_adaptive || mesh_point_bool) {
1338
1339 // Update temporaries
1340 step += 1;
1341 h_step = h_tmp_step;
1342
1343 // Update outputs
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);
1347
1348 // Callback after the step is completed
1349 if (this->m_step_callback) {
1350 this->m_step_callback(step, x_new_step, t_step);
1351 }
1352
1353 // Check if the current step is the last one
1354 if (std::abs(t_step - t_mesh(last)) < SQRT_EPSILON) {break;}
1355
1356 // Update the previous step
1357 x_old_step = x_new_step;
1358 }
1359 }
1360 return true;
1361 }
1362
1374 bool adaptive_solve(VectorX const & t_mesh, VectorN const & ics, Solution<Real, N, M> & sol) const
1375 {
1376 using Eigen::all;
1377 using Eigen::last;
1378
1379 #define CMD "Sandals::RungeKutta::adaptive_solve(...): "
1380
1381 // Check if the adaptive method is enabled and the method is embedded
1382 if (!this->is_embedded()) {
1383 SANDALS_WARNING(CMD "the method is not embedded, using solve(...) method.");
1384 return this->solve(t_mesh, ics, sol);
1385 } else if (!this->m_adaptive) {
1386 SANDALS_WARNING(CMD "adaptive method is disabled, using solve(...) method.");
1387 return this->solve(t_mesh, ics, sol);
1388 }
1389
1390 // Instantiate output
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};
1393 if (this->m_tableau.is_embedded) {
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);
1396 } else {
1397 sol.resize(t_mesh.size());
1398 }
1399
1400 // Store initial conditions
1401 sol.t(0) = t_mesh(0);
1402 sol.x.col(0) = ics;
1403 sol.h.col(0) = this->m_system->h(ics, t_mesh(0));
1404
1405 // Callback on initial conditions
1406 if (this->m_step_callback) {this->m_step_callback(0, ics, t_mesh(0));}
1407
1408 // Instantiate temporary variables
1409 Integer step{0};
1410 VectorN x_old_step(ics), x_new_step(ics);
1411
1412 while (true) {
1413 // Integrate system
1414 this->advance(x_old_step, t_step, h_step, x_new_step, h_new_step);
1415
1416 // Update the current step
1417 t_step += h_step;
1418
1419 // Saturate the suggested timestep
1420 if (this->m_adaptive && this->m_tableau.is_embedded) {
1421 h_step = std::max(std::min(h_new_step, h_max), h_min);
1422 }
1423
1424 SANDALS_ASSERT(step < sol.size(), CMD "safety length exceeded.");
1425
1426 // Update temporaries
1427 step += 1;
1428
1429 // Update outputs
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);
1433
1434 // Callback after the step is completed
1435 if (this->m_step_callback) {
1436 this->m_step_callback(step, x_new_step, t_step);
1437 }
1438
1439 // Check if the current step is the last one
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;}
1442
1443 // Update the previous step
1444 x_old_step = x_new_step;
1445 }
1446
1447 // Resize the output
1449
1450 return true;
1451
1452 #undef CMD
1453 }
1454
1464 bool project(VectorN const & x, Real const t, VectorN & x_projected) const
1465 {
1466 #define CMD "Sandals::RungeKutta::project(...): "
1467
1468 // Check if there are any constraints
1469 x_projected = x;
1470 if (M > Integer(0)) {
1471
1472 VectorM h;
1473 MatrixM Jh_x;
1474 VectorP b, x_step;
1475 MatrixP A;
1476 A.setZero();
1477 A.template block<N, N>(0, 0) = MatrixN::Identity();
1478 for (Integer k{0}; k < this->m_max_projection_iterations; ++k) {
1479
1480 /* Standard projection method
1481 [A] {x} = {b}
1482 / I Jh_x^T \ / dx \ / x_t - x_k \
1483 | | | | = | |
1484 \ Jh_x 0 / \ lambda / \ -h /
1485 */
1486
1487 // Evaluate the invariants vector and its Jacobian
1488 h = this->m_system->h(x_projected, t);
1489 Jh_x = this->m_system->Jh_x(x_projected, t);
1490
1491 // Check if the solution is found
1492 if (h.norm() < this->m_projection_tolerance) {return true;}
1493
1494 // Build the linear system
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;
1499
1500 // Compute the solution of the linear system
1501 this->m_lu.compute(A);
1502 SANDALS_ASSERT(this->m_lu.rank() == N+M, CMD "singular Jacobian detected.");
1503 x_step = this->m_lu.solve(b);
1504
1505 // Check if the step is too small
1506 if (x_step.norm() < this->m_projection_tolerance * this->m_projection_tolerance) {return false;}
1507
1508 // Update the solution
1509 x_projected.noalias() += x_step(Eigen::seqN(0, N));
1510 }
1511 if (this->m_verbose) {SANDALS_WARNING(CMD "maximum number of iterations reached.");}
1512 return false;
1513 } else {
1514 return true;
1515 }
1516
1517 #undef CMD
1518 }
1519
1531 bool project_ics(VectorN const & x, Real const t, std::vector<Integer> const & projected_equations,
1532 std::vector<Integer> const & projected_invariants, VectorN & x_projected) const
1533 {
1534 #define CMD "Sandals::RungeKutta::project_ics(...): "
1535
1536 Integer X{static_cast<Integer>(projected_equations.size())};
1537 Integer H{static_cast<Integer>(projected_invariants.size())};
1538
1539 // Check if there are any constraints
1540 x_projected = x;
1541 if (H > Integer(0)) {
1542
1543 VectorM h;
1544 MatrixM Jh_x;
1545 VectorX b(X+H), x_step(X+H);
1546 MatrixX A(X+H, X+H);
1547 A.setZero();
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) {
1551
1552 /* Standard projection method
1553 [A] {x} = {b}
1554 / I Jh_x^T \ / dx \ / x_t - x_k \
1555 | | | | = | |
1556 \ Jh_x 0 / \ lambda / \ -h /
1557 */
1558
1559 // Evaluate the invariants vector and its Jacobian
1560 h = this->m_system->h(x_projected, t);
1561 Jh_x = this->m_system->Jh_x(x_projected, t);
1562
1563 // Select only the projected invariants
1564 h = h(projected_invariants);
1565 Jh_x = Jh_x(projected_invariants, projected_equations);
1566
1567 // Check if the solution is found
1568 if (h.norm() < this->m_projection_tolerance) {return true;}
1569
1570 // Build the linear system
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);
1574 b.tail(H) = -h;
1575
1576 // Compute the solution of the linear system
1577 lu.compute(A);
1578 SANDALS_ASSERT(lu.rank() == X+H, CMD "singular Jacobian detected.");
1579 x_step = this->m_lu.solve(b);
1580
1581 // Check if the step is too small
1582 if (x_step.norm() < this->m_projection_tolerance * this->m_projection_tolerance) {return false;}
1583
1584 // Update the solution
1585 x_projected(projected_equations).noalias() += x_step;
1586 }
1587 if (this->m_verbose) {SANDALS_WARNING(CMD "maximum number of iterations reached.");}
1588 return false;
1589 } else {
1590 return true;
1591 }
1592
1593 #undef CMD
1594 }
1595
1603 Real estimate_order(std::vector<VectorX> const & t_mesh, VectorN const & ics, std::function<MatrixX(VectorX)> & sol) const
1604 {
1605 using Eigen::last;
1606
1607 #define CMD "Sandals::RungeKutta::estimate_order(...): "
1608
1609 SANDALS_ASSERT(t_mesh.size() > Integer(1), CMD "expected at least two time meshes.");
1610
1611 for (Integer i{0}; i < static_cast<Integer>(t_mesh.size()); ++i) {
1612
1613 // Check if the time mesh has the same initial and final time
1614 SANDALS_ASSERT((t_mesh[0](0) - t_mesh[i](0)) < SQRT_EPSILON,
1615 CMD "expected the same initial time.");
1616 SANDALS_ASSERT((t_mesh[0](last) - t_mesh[i](last)) < SQRT_EPSILON,
1617 CMD "expected the same final time.");
1618
1619 // Check if the mesh step is fixed all over the time mesh
1620 for (Integer j{1}; j < static_cast<Integer>(t_mesh[i].size()); ++j) {
1621 SANDALS_ASSERT((t_mesh[i](j) - t_mesh[i](j-1)) - (t_mesh[i](1) - t_mesh[i](0)) < SQRT_EPSILON,
1622 CMD "expected a fixed step.");
1623 }
1624 }
1625
1626 // Solve the system for each time scale
1627 Solution<Real, N, M> sol_num;
1628 MatrixX sol_ana;
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) {
1631 SANDALS_ASSERT(this->solve(t_mesh[i], ics, sol_num), CMD "failed to solve the system for " <<
1632 "the" << i << "-th time mesh.");
1633 sol_ana = sol(sol_num.t);
1634 SANDALS_ASSERT(sol_ana.rows() == sol_num.x.rows(),
1635 CMD "expected the same number of states in analytical solution.");
1636 SANDALS_ASSERT(sol_ana.cols() == sol_num.x.cols(),
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();
1640 }
1641
1642 // Compute the order of the method thorugh least squares
1643 VectorX A(h_vec.array().log());
1644 VectorX b(e_vec.array().log());
1645 return ((A.transpose() * A).ldlt().solve(A.transpose() * b))(0);
1646
1647 #undef CMD
1648 }
1649
1650 }; // class RungeKutta
1651
1652} // namespace Sandals
1653
1654#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: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