update mini-odeint with tests
This commit is contained in:
9
mini-odeint/.gitignore
vendored
Normal file
9
mini-odeint/.gitignore
vendored
Normal file
@@ -0,0 +1,9 @@
|
||||
.DS_Store
|
||||
.idea
|
||||
*.log
|
||||
tmp/
|
||||
|
||||
.cache/
|
||||
compile_commands.json
|
||||
/build/
|
||||
*.o
|
||||
28
mini-odeint/CMakeLists.txt
Normal file
28
mini-odeint/CMakeLists.txt
Normal file
@@ -0,0 +1,28 @@
|
||||
cmake_minimum_required(VERSION 3.20)
|
||||
|
||||
project(mini-odeint)
|
||||
|
||||
option(FORCE_FETCH_CATCH2 "Force fetching Catch2" OFF)
|
||||
|
||||
set(CMAKE_CXX_STANDARD 20)
|
||||
set(CMAKE_CXX_STANDARD_REQUIRED ON)
|
||||
|
||||
find_package(Catch2 3 QUIET)
|
||||
|
||||
if(NOT TARGET Catch2::Catch2WithMain OR FORCE_FETCH_CATCH2)
|
||||
include(FetchContent)
|
||||
FetchContent_Declare(
|
||||
Catch2
|
||||
GIT_REPOSITORY https://github.com/catchorg/Catch2.git
|
||||
GIT_TAG v3.4.0)
|
||||
FetchContent_MakeAvailable(Catch2)
|
||||
list(APPEND CMAKE_MODULE_PATH ${catch2_SOURCE_DIR}/extras)
|
||||
endif()
|
||||
|
||||
include(CTest)
|
||||
include(Catch)
|
||||
|
||||
add_executable(tests mini-odeint-tests.cpp mini-odeint.hpp)
|
||||
target_link_libraries(tests PRIVATE Catch2::Catch2WithMain)
|
||||
|
||||
catch_discover_tests(tests)
|
||||
126
mini-odeint/mini-odeint-tests.cpp
Normal file
126
mini-odeint/mini-odeint-tests.cpp
Normal file
@@ -0,0 +1,126 @@
|
||||
#include <catch2/catch_test_macros.hpp>
|
||||
#include <catch2/matchers/catch_matchers_floating_point.hpp>
|
||||
#include <catch2/matchers/catch_matchers_range_equals.hpp>
|
||||
|
||||
#include "mini-odeint.hpp"
|
||||
|
||||
#include <ranges>
|
||||
|
||||
unsigned int Factorial(unsigned int number) {
|
||||
return number <= 1 ? number : Factorial(number - 1) * number;
|
||||
}
|
||||
|
||||
TEST_CASE("Factorials are computed", "[factorial]") {
|
||||
REQUIRE(Factorial(1) == 1);
|
||||
REQUIRE(Factorial(2) == 2);
|
||||
REQUIRE(Factorial(3) == 6);
|
||||
REQUIRE(Factorial(10) == 3628800);
|
||||
}
|
||||
|
||||
TEST_CASE("Mini-odeint works", "[mini-odeint]") {
|
||||
using namespace mini_odeint;
|
||||
|
||||
SECTION("float") {
|
||||
std::vector<float> times;
|
||||
times.reserve(1001);
|
||||
for (int i = 0; i <= 1000; ++i) {
|
||||
times.push_back(i / 1000.0f);
|
||||
}
|
||||
|
||||
SECTION("scalar") {
|
||||
std::vector<float> ys(times.size());
|
||||
|
||||
explicitRungeKutta(std::span(ys), std::span<const float>(times), 1.0f,
|
||||
1e-6f, [](auto y, auto t) { return -y; });
|
||||
|
||||
REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
|
||||
return std::abs(a - std::exp(-b)) < 1e-6;
|
||||
}));
|
||||
}
|
||||
|
||||
SECTION("2 vector") {
|
||||
std::vector<std::array<float, 2>> ys(times.size());
|
||||
|
||||
explicitRungeKutta(
|
||||
std::span(ys), std::span<const float>(times), std::array{1.0f, 2.0f},
|
||||
1e-6f, [](auto y, auto t) { return std::array{-y[0], -y[1]}; });
|
||||
|
||||
REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
|
||||
const auto s = std::exp(-b);
|
||||
|
||||
return func::inf_norm(std::array{a[0] - s, a[1] - 2 * s}) <
|
||||
1e-6;
|
||||
}));
|
||||
}
|
||||
|
||||
SECTION("4 vector") {
|
||||
std::vector<std::array<float, 4>> ys(times.size());
|
||||
|
||||
explicitRungeKutta(std::span(ys), std::span<const float>(times),
|
||||
std::array{1.0f, 2.0f, 3.0f, 4.0f}, 1e-6f,
|
||||
[](auto y, auto t) {
|
||||
return std::array{-y[0], -y[1], -y[2], -y[3]};
|
||||
});
|
||||
|
||||
REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
|
||||
const auto s = std::exp(-b);
|
||||
|
||||
return func::inf_norm(std::array{a[0] - s, a[1] - 2 * s,
|
||||
a[2] - 3 * s,
|
||||
a[3] - 4 * s}) < 1e-6;
|
||||
}));
|
||||
}
|
||||
}
|
||||
|
||||
SECTION("float") {
|
||||
std::vector<double> times;
|
||||
times.reserve(1001);
|
||||
for (int i = 0; i <= 1000; ++i) {
|
||||
times.push_back(i / 1000.0);
|
||||
}
|
||||
|
||||
SECTION("scalar") {
|
||||
std::vector<double> ys(times.size());
|
||||
|
||||
explicitRungeKutta(std::span(ys), std::span<const double>(times), 1.0,
|
||||
1e-6, [](auto y, auto t) { return -y; });
|
||||
|
||||
REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
|
||||
return std::abs(a - std::exp(-b)) < 1e-6;
|
||||
}));
|
||||
}
|
||||
|
||||
SECTION("2 vector") {
|
||||
std::vector<std::array<double, 2>> ys(times.size());
|
||||
|
||||
explicitRungeKutta(
|
||||
std::span(ys), std::span<const double>(times), std::array{1.0, 2.0},
|
||||
1e-6, [](auto y, auto t) { return std::array{-y[0], -y[1]}; });
|
||||
|
||||
REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
|
||||
const auto s = std::exp(-b);
|
||||
|
||||
return func::inf_norm(std::array{a[0] - s, a[1] - 2 * s}) <
|
||||
1e-6;
|
||||
}));
|
||||
}
|
||||
|
||||
SECTION("4 vector") {
|
||||
std::vector<std::array<double, 4>> ys(times.size());
|
||||
|
||||
explicitRungeKutta(std::span(ys), std::span<const double>(times),
|
||||
std::array{1.0, 2.0, 3.0, 4.0}, 1e-6,
|
||||
[](auto y, auto t) {
|
||||
return std::array{-y[0], -y[1], -y[2], -y[3]};
|
||||
});
|
||||
|
||||
REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
|
||||
const auto s = std::exp(-b);
|
||||
|
||||
return func::inf_norm(std::array{a[0] - s, a[1] - 2 * s,
|
||||
a[2] - 3 * s,
|
||||
a[3] - 4 * s}) < 1e-6;
|
||||
}));
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -4,8 +4,6 @@
|
||||
#include <array>
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
#include <iterator>
|
||||
#include <ranges>
|
||||
#include <span>
|
||||
|
||||
namespace mini_odeint {
|
||||
@@ -64,6 +62,127 @@ template <typename T> struct DormandPrince {
|
||||
90730570.0 / 29380423.0, -8293050.0 / 29380423.0}}};
|
||||
};
|
||||
|
||||
template <typename T> struct scalar_type {
|
||||
using type = T;
|
||||
};
|
||||
|
||||
template <typename T, std::size_t N> struct scalar_type<T[N]> {
|
||||
using type = T;
|
||||
};
|
||||
|
||||
template <typename T, std::size_t N> struct scalar_type<std::array<T, N>> {
|
||||
using type = T;
|
||||
};
|
||||
|
||||
template <typename T> using scalar_type_t = typename scalar_type<T>::type;
|
||||
|
||||
template <typename T> struct is_std_array : std::false_type {};
|
||||
template <typename T, std::size_t N>
|
||||
struct is_std_array<std::array<T, N>> : std::true_type {};
|
||||
template <typename T>
|
||||
inline constexpr bool is_std_array_v = is_std_array<T>::value;
|
||||
|
||||
template <typename T>
|
||||
requires std::is_trivially_copyable_v<T>
|
||||
class OdeVector {
|
||||
|
||||
using Scalar = scalar_type_t<T>;
|
||||
|
||||
T value;
|
||||
|
||||
public:
|
||||
OdeVector() = default;
|
||||
explicit OdeVector(T value) : value(std::move(value)) {}
|
||||
|
||||
friend constexpr OdeVector operator*(OdeVector lhs, const Scalar &rhs) {
|
||||
return lhs *= rhs;
|
||||
}
|
||||
|
||||
friend constexpr OdeVector operator*(const Scalar &lhs, OdeVector rhs) {
|
||||
return rhs *= lhs;
|
||||
}
|
||||
|
||||
friend constexpr OdeVector operator+(OdeVector lhs, const OdeVector &rhs) {
|
||||
return lhs += rhs;
|
||||
}
|
||||
|
||||
friend constexpr OdeVector operator+(const T &lhs, OdeVector rhs) {
|
||||
return rhs += OdeVector{lhs};
|
||||
}
|
||||
|
||||
friend constexpr OdeVector operator+(OdeVector lhs, const T &rhs) {
|
||||
return lhs += OdeVector{rhs};
|
||||
}
|
||||
|
||||
const OdeVector &operator*=(const Scalar &rhs) {
|
||||
if constexpr (std::is_bounded_array_v<T>) {
|
||||
for (std::size_t i = 0; i < std::extent_v<T>; ++i) {
|
||||
value[i] *= rhs;
|
||||
}
|
||||
} else if constexpr (is_std_array_v<T>) {
|
||||
for (std::size_t i = 0; i < std::tuple_size_v<T>; ++i) {
|
||||
value[i] *= rhs;
|
||||
}
|
||||
} else {
|
||||
value *= rhs;
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
constexpr OdeVector &operator+=(const OdeVector &rhs) {
|
||||
if constexpr (std::is_bounded_array_v<T>) {
|
||||
for (std::size_t i = 0; i < std::extent_v<T>; ++i) {
|
||||
value[i] += rhs.value[i];
|
||||
}
|
||||
} else if constexpr (is_std_array_v<T>) {
|
||||
for (std::size_t i = 0; i < std::tuple_size_v<T>; ++i) {
|
||||
value[i] += rhs.value[i];
|
||||
}
|
||||
} else {
|
||||
value += rhs.value;
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
constexpr OdeVector &operator=(const T &rhs) {
|
||||
value = rhs;
|
||||
return *this;
|
||||
}
|
||||
|
||||
explicit constexpr operator const T &() const { return value; }
|
||||
};
|
||||
|
||||
template <typename E> struct Vec2 {
|
||||
using value_type = E;
|
||||
|
||||
E x, y;
|
||||
|
||||
constexpr Vec2() = default;
|
||||
constexpr Vec2(E x, E y) : x(x), y(y) {}
|
||||
|
||||
constexpr explicit Vec2(std::array<E, 3> v) : x(v[0]), y(v[1]) {}
|
||||
|
||||
friend constexpr Vec2 operator+(Vec2 lhs, const Vec2 &rhs) {
|
||||
return lhs += rhs;
|
||||
}
|
||||
friend constexpr Vec2 operator*(Vec2 lhs, const value_type &rhs) {
|
||||
return lhs *= rhs;
|
||||
}
|
||||
friend constexpr Vec2 operator*(const value_type &lhs, Vec2 rhs) {
|
||||
return rhs *= lhs;
|
||||
}
|
||||
constexpr Vec2 &operator+=(const Vec2 &rhs) {
|
||||
x += rhs.x;
|
||||
y += rhs.y;
|
||||
return *this;
|
||||
}
|
||||
constexpr Vec2 &operator*=(const value_type &rhs) {
|
||||
x *= rhs;
|
||||
y *= rhs;
|
||||
return *this;
|
||||
}
|
||||
};
|
||||
|
||||
template <typename E> struct Vec3 {
|
||||
using value_type = E;
|
||||
|
||||
@@ -74,14 +193,14 @@ template <typename E> struct Vec3 {
|
||||
|
||||
constexpr explicit Vec3(std::array<E, 3> v) : x(v[0]), y(v[1]), z(v[2]) {}
|
||||
|
||||
friend constexpr Vec3 operator+(const Vec3 &lhs, const Vec3 &rhs) {
|
||||
return Vec3{lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z};
|
||||
friend constexpr Vec3 operator+(Vec3 lhs, const Vec3 &rhs) {
|
||||
return lhs += rhs;
|
||||
}
|
||||
friend constexpr Vec3 operator*(const Vec3 &lhs, const value_type &rhs) {
|
||||
return Vec3{lhs.x * rhs, lhs.y * rhs, lhs.z * rhs};
|
||||
friend constexpr Vec3 operator*(Vec3 lhs, const value_type &rhs) {
|
||||
return lhs *= rhs;
|
||||
}
|
||||
friend constexpr Vec3 operator*(const value_type &lhs, const Vec3 &rhs) {
|
||||
return Vec3{lhs * rhs.x, lhs * rhs.y, lhs * rhs.z};
|
||||
friend constexpr Vec3 operator*(const value_type &lhs, Vec3 rhs) {
|
||||
return rhs *= lhs;
|
||||
}
|
||||
constexpr Vec3 &operator+=(const Vec3 &rhs) {
|
||||
x += rhs.x;
|
||||
@@ -89,26 +208,41 @@ template <typename E> struct Vec3 {
|
||||
z += rhs.z;
|
||||
return *this;
|
||||
}
|
||||
constexpr Vec3 &operator*=(const value_type &rhs) {
|
||||
x *= rhs;
|
||||
y *= rhs;
|
||||
z *= rhs;
|
||||
return *this;
|
||||
}
|
||||
};
|
||||
|
||||
namespace alg {
|
||||
template <typename T> inline auto inf_norm(const T &v) {
|
||||
std::ranges::max_element(v, {}, [](auto n) { return std::abs(n); });
|
||||
namespace func {
|
||||
|
||||
inline auto inf_norm(std::ranges::range auto v) {
|
||||
return *std::ranges::max_element(v, {}, [](auto n) { return std::abs(n); });
|
||||
}
|
||||
|
||||
inline float inf_norm(float v) { return std::abs(v); }
|
||||
inline std::floating_point auto inf_norm(std::floating_point auto v) {
|
||||
return std::abs(v);
|
||||
}
|
||||
|
||||
inline double inf_norm(double v) { return std::abs(v); }
|
||||
template <typename T> inline auto inf_norm(const OdeVector<T> &v) {
|
||||
return inf_norm(static_cast<const T &>(v));
|
||||
}
|
||||
|
||||
template <typename E> inline E inf_norm(const Vec2<E> &v) {
|
||||
return std::max(std::abs(v.x), std::abs(v.y));
|
||||
}
|
||||
|
||||
template <typename E> inline E inf_norm(const Vec3<E> &v) {
|
||||
return std::max({std::abs(v.x), std::abs(v.y), std::abs(v.z)});
|
||||
}
|
||||
|
||||
} // namespace alg
|
||||
} // namespace func
|
||||
|
||||
template <typename Vector, typename Scalar = std::iter_value_t<Vector>,
|
||||
template <typename Vector, typename Scalar = scalar_type_t<Vector>,
|
||||
typename Tableau = DormandPrince<Scalar>>
|
||||
requires std::same_as<Scalar, std::iter_value_t<Tableau>>
|
||||
requires std::same_as<Scalar, typename Tableau::value_type>
|
||||
|
||||
inline std::size_t explicitRungeKutta(std::span<Vector> ys,
|
||||
std::span<Scalar const> ts, Vector y0,
|
||||
@@ -126,14 +260,14 @@ inline std::size_t explicitRungeKutta(std::span<Vector> ys,
|
||||
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");
|
||||
static_assert(c.back() == 1.0, "last c value must be 1.0");
|
||||
|
||||
auto y_hat_n = y0;
|
||||
OdeVector<Vector> y_hat_n{y0};
|
||||
ys[0] = y0;
|
||||
|
||||
std::size_t it = 1;
|
||||
|
||||
std::array<Vector, stages> k;
|
||||
std::array<OdeVector<Vector>, stages> k;
|
||||
|
||||
const auto N = ts.size();
|
||||
if (!N) {
|
||||
@@ -143,8 +277,6 @@ inline std::size_t explicitRungeKutta(std::span<Vector> ys,
|
||||
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]) {
|
||||
@@ -154,16 +286,17 @@ inline std::size_t explicitRungeKutta(std::span<Vector> ys,
|
||||
const auto last_k_store = k[stages - 1];
|
||||
k[0] = k[stages - 1];
|
||||
for (std::size_t i = 1; i < stages; ++i) {
|
||||
Vector sum_ak{};
|
||||
OdeVector<Vector> 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);
|
||||
k[i] =
|
||||
dydx(static_cast<Vector>(y_hat_n + h_n * sum_ak), t_n + c[i] * h_n);
|
||||
}
|
||||
|
||||
// calculate final value and error
|
||||
Vector error{};
|
||||
Vector sum_bk{};
|
||||
OdeVector<Vector> error{};
|
||||
OdeVector<Vector> 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];
|
||||
@@ -171,13 +304,13 @@ inline std::size_t explicitRungeKutta(std::span<Vector> ys,
|
||||
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 = alg::inf_norm(h_n * error);
|
||||
const auto E_hp1 = func::inf_norm(h_n * error);
|
||||
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;
|
||||
Vector Phi{};
|
||||
OdeVector<Vector> Phi{};
|
||||
for (std::size_t i = 0; i < stages; ++i) {
|
||||
auto term = sigma;
|
||||
auto b_i = term * p[i][0];
|
||||
@@ -187,7 +320,7 @@ inline std::size_t explicitRungeKutta(std::span<Vector> ys,
|
||||
}
|
||||
Phi += b_i * k[i];
|
||||
}
|
||||
ys[it] = y_hat_n + h_n * Phi;
|
||||
ys[it] = static_cast<Vector>(y_hat_n + h_n * Phi);
|
||||
++it;
|
||||
}
|
||||
|
||||
@@ -195,7 +328,6 @@ inline std::size_t explicitRungeKutta(std::span<Vector> ys,
|
||||
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;
|
||||
@@ -212,28 +344,4 @@ inline std::size_t explicitRungeKutta(std::span<Vector> ys,
|
||||
|
||||
} // namespace mini_odeint
|
||||
|
||||
#include <vector>
|
||||
|
||||
#include <iostream>
|
||||
|
||||
int main() {
|
||||
using namespace mini_odeint;
|
||||
// make vector of floats from 0.0 to 1.0 by 0.001
|
||||
std::vector<float> times;
|
||||
times.reserve(1001);
|
||||
for (int i = 0; i < 1000; ++i) {
|
||||
times.push_back(i / 1000.0);
|
||||
}
|
||||
|
||||
std::vector<Vec3<float>> ys(times.size());
|
||||
|
||||
explicitRungeKutta(
|
||||
std::span(ys), std::span<const float>(times), Vec3<float>{1.0, 1.0, 1.0},
|
||||
float(1e-6), [](auto y, auto t) { return Vec3<float>{-y.x, 0.0, 0.0}; });
|
||||
|
||||
for (const auto &v : ys) {
|
||||
std::cout << v.x << '\n';
|
||||
}
|
||||
}
|
||||
|
||||
#endif // MINI_ODEINT_H_
|
||||
Reference in New Issue
Block a user