From 02b3a309e1502a820ace8857b68ceb4800cec2d8 Mon Sep 17 00:00:00 2001 From: "John K. Luebs" Date: Thu, 7 Nov 2024 22:42:55 -0600 Subject: [PATCH] Initial --- mini-odeint/.clangd | 2 + mini-odeint/mini-odeint.cpp | 288 ++++++++++++++++++++++++++++++++++++ 2 files changed, 290 insertions(+) create mode 100644 mini-odeint/.clangd create mode 100644 mini-odeint/mini-odeint.cpp diff --git a/mini-odeint/.clangd b/mini-odeint/.clangd new file mode 100644 index 0000000..359a391 --- /dev/null +++ b/mini-odeint/.clangd @@ -0,0 +1,2 @@ +CompileFlags: + Add: [-std=c++20] diff --git a/mini-odeint/mini-odeint.cpp b/mini-odeint/mini-odeint.cpp new file mode 100644 index 0000000..16c5742 --- /dev/null +++ b/mini-odeint/mini-odeint.cpp @@ -0,0 +1,288 @@ +#ifndef MINI_ODEINT_H_ +#define MINI_ODEINT_H_ + +#include +#include +#include +#include +#include + +namespace mini_odeint { + +template struct DormandPrince { + using value_type = T; + static constexpr std::size_t stages = 7; + static constexpr std::size_t order = 5; + static constexpr std::size_t estimator_order = 4; + static constexpr std::size_t dense_order = 5; + + static constexpr std::array c{ + 0.0, 1.0 / 5.0, 3.0 / 10.0, 4.0 / 5.0, 8.0 / 9.0, 1.0, 1.0}; + + static constexpr std::array, stages> a{ + {{0.0, 0.0, 0.0, 0.0, 0.0}, + {1.0 / 5.0, 0.0, 0.0, 0.0, 0.0}, + {3.0 / 40.0, 9.0 / 40.0, 0.0, 0.0, 0.0}, + {44.0 / 45.0, -56.0 / 15.0, 32.0 / 9.0, 0.0, 0.0}, + {19372.0 / 6561.0, -25360.0 / 2187.0, 64448.0 / 6561.0, -212.0 / 729.0, + 0.0}, + {9017.0 / 3168.0, -355.0 / 33.0, 46732.0 / 5247.0, 49.0 / 176.0, + -5103.0 / 18656.0}, + {35.0 / 384.0, 0.0, 500.0 / 1113.0, 125.0 / 192.0, -2187.0 / 6784.0, + 11.0 / 84.0}}}; + + static constexpr std::array b_hat{ + 35.0 / 384.0, 0.0, 500.0 / 1113.0, 125.0 / 192.0, -2187.0 / 6784.0, + 11.0 / 84.0, 0.0}; + + static constexpr std::array b{5179.0 / 57600.0, 0.0, + 7571.0 / 16695.0, 393.0 / 640.0, + -92097.0 / 339200.0, 187.0 / 2100.0, + 1.0 / 40.0}; + + static constexpr std::array, stages> p{ + + {{1.0, -32272833064.0 / 11282082432.0, 34969693132.0 / 11282082432.0, + -13107642775.0 / 11282082432.0, 157015080.0 / 11282082432.0}, + {0.0, 0.0, 0.0, 0.0, 0.0}, + {0.0, 1323431896.0 * 100.0 / 32700410799.0, + -2074956840.0 * 100.0 / 32700410799.0, + 914128567.0 * 100.0 / 32700410799.0, + -15701508.0 * 100.0 / 32700410799.0}, + {0.0, -889289856.0 * 25.0 / 5641041216.0, + 2460397220.0 * 25.0 / 5641041216.0, -1518414297.0 * 25.0 / 5641041216.0, + 94209048.0 * 25.0 / 5641041216.0}, + {0.0, 259006536.0 * 2187.0 / 199316789632.0, + -687873124.0 * 2187.0 / 199316789632.0, + 451824525.0 * 2187.0 / 199316789632.0, + -52338360.0 * 2187.0 / 199316789632.0}, + {0.0, -361440756.0 * 11.0 / 2467955532.0, + 946554244.0 * 11.0 / 2467955532.0, -661884105.0 * 11.0 / 2467955532.0, + 106151040.0 * 11.0 / 2467955532.0}, + {0.0, 44764047.0 / 29380423.0, -127201567 / 29380423.0, + 90730570.0 / 29380423.0, -8293050.0 / 29380423.0}}}; +}; + +template struct OdeVector { + using value_type = E; + + E value; + + constexpr OdeVector() = default; + constexpr explicit OdeVector(E value) : value(std::move(value)) {} + + constexpr operator E() const { return value; } + + friend constexpr OdeVector operator+(const OdeVector &lhs, + const OdeVector &rhs) { + return OdeVector{lhs.value + rhs.value}; + } + friend constexpr OdeVector operator*(const OdeVector &lhs, + const value_type &rhs) { + return OdeVector{lhs.value * rhs}; + } + friend constexpr OdeVector operator*(const value_type &lhs, + const OdeVector &rhs) { + return OdeVector{lhs * rhs.value}; + } + constexpr OdeVector &operator+=(const OdeVector &rhs) { + value += rhs.value; + return *this; + } + constexpr OdeVector &operator+=(const value_type &rhs) { + value += rhs; + return *this; + } + value_type inf_norm() const { return std::abs(value); } +}; + +template struct OdeVector> { + using array_type = std::array; + using value_type = T; + std::array value; + + constexpr OdeVector() = default; + constexpr explicit OdeVector(std::array value) + : value(std::move(value)) {} + + operator std::array() const { return value; } + + friend constexpr OdeVector operator+(const OdeVector &lhs, + const OdeVector &rhs) { + OdeVector result{lhs}; + for (std::size_t i = 0; i < N; ++i) { + result.value[i] += rhs.value[i]; + } + return result; + } + friend constexpr OdeVector operator*(const OdeVector &lhs, + const value_type &rhs) { + OdeVector result{lhs}; + for (std::size_t i = 0; i < N; ++i) { + result.value[i] *= rhs; + } + return result; + } + friend constexpr OdeVector operator*(const value_type &lhs, + const OdeVector &rhs) { + OdeVector result{rhs}; + for (std::size_t i = 0; i < N; ++i) { + result.value[i] *= lhs; + } + return result; + } + constexpr OdeVector &operator+=(const OdeVector &rhs) { + for (std::size_t i = 0; i < N; ++i) { + value[i] += rhs.value[i]; + } + return *this; + } + + constexpr OdeVector &operator+=(const value_type &rhs) { + for (std::size_t i = 0; i < N; ++i) { + value[i] += rhs; + } + return *this; + } + + value_type inf_norm() const { + auto max = std::abs(value[0]); + for (std::size_t i = 1; i < N; ++i) { + max = std::max(max, std::abs(value[i])); + } + return max; + } +}; + +template , + typename Tableau = DormandPrince> + requires std::same_as> + +inline std::size_t explicitRungeKutta(std::span ys, + std::span ts, Vector y0, + Scalar tol, auto &&dydx) + requires requires(decltype(dydx) f) { + { f(Vector{}, Scalar{}) } -> std::same_as; + } +{ + const auto stages = Tableau::stages; + const auto dense_order = Tableau::dense_order; + const auto order = Tableau::order; + const auto &a = Tableau::a; + const auto &p = Tableau::p; + const auto &c = Tableau::c; + const auto &b = Tableau::b; + const auto &b_hat = Tableau::b_hat; + + static_assert(Tableau::c.back() == 1.0, "last c value must be 1.0"); + + auto y_hat_n = y0; + ys[0] = y0; + + std::size_t it = 1; + + std::array k; + + const auto N = ts.size(); + if (!N) { + return 0; + } + + auto t_n = ts[0]; + auto h_n = ts[N - 1] - t_n; + + int step_count = 0; + + k[stages - 1] = dydx(y0, t_n); + + while (t_n < ts[N - 1]) { + auto step_rejected = true; + while (step_rejected) { + // reuse last k (we have asserted that the last c value is 1.0) + const auto last_k_store = k[stages - 1]; + k[0] = k[stages - 1]; + for (std::size_t i = 1; i < stages; ++i) { + OdeVector sum_ak{}; + for (std::size_t j = 0; j < i; ++j) { + sum_ak += a[i][j] * k[j]; + } + k[i] = dydx(y_hat_n + h_n * sum_ak, t_n + c[i] * h_n); + } + + // calculate final value and error + OdeVector error{}; + OdeVector sum_bk{}; + for (std::size_t i = 0; i < stages; ++i) { + sum_bk += b_hat[i] * k[i]; + error += (b_hat[i] - b[i]) * k[i]; + } + const auto y_hat_np1 = y_hat_n + h_n * sum_bk; + + // check if step is successful, ie error is within tolerance + const auto E_hp1 = (h_n * error).inf_norm(); + if (E_hp1 < tol) { + // if moved over any requested times then interpolate their values + const auto t_np1 = t_n + h_n; + while (it < N && t_np1 >= ts[it]) { + const auto sigma = (ts[it] - t_n) / h_n; + OdeVector Phi{}; + for (std::size_t i = 0; i < stages; ++i) { + auto term = sigma; + auto b_i = term * p[i][0]; + for (std::size_t j = 1; j < dense_order; ++j) { + term *= sigma; + b_i += term * p[i][j]; + } + Phi += b_i * k[i]; + } + ys[it] = y_hat_n + h_n * Phi; + ++it; + } + + // move to next step + step_rejected = false; + y_hat_n = y_hat_np1; + t_n = t_np1; + ++step_count; + } else { + // failed step, reset last k back to stored value + k[stages - 1] = last_k_store; + } + + // adapt step size + h_n *= 0.9 * std::pow(tol / E_hp1, 1.0 / (order + 1.0)); + } + } + + assert(it == N); + return it; +} + +} // namespace mini_odeint + +#include + +#include + +int main() { + using namespace mini_odeint; + // make vector of doubles from 0.0 to 1.0 by 0.001 + std::vector times; + times.reserve(1001); + for (int i = 0; i < 1000; ++i) { + times.push_back(i / 1000.0); + } + + std::vector> ys(times.size()); + + explicitRungeKutta( + std::span(ys), std::span(times), + std::array{1.0, 2.0}, float(1e-6), + [](auto y, auto t) { return std::array{-y[0], -y[1]}; }); + + for (const auto &v : ys) { + std::cout << v[0] << v[1] << '\n'; + } +} + +#endif // MINI_ODEINT_H_