Instantly share code, notes, and snippets.

Embed
What would you like to do?
#include <iostream>
#include <fstream>
#include <array>
#include <vector>
#include "cranberries/rkf45.hpp"
int main()
{
using state_type = std::array<long double, 3>;
std::vector<long double> timelog;
std::vector<state_type> statelog;
state_type state0 = { { 10, 1, 1 } };
bool exit_status = cranberries::make_rkf45(
// Lorenz attractor
[&, p = 10.L, r = 28.L, b = 8.L / 3](const long double t, const state_type& x) {
state_type dxdt{};
// x=x[0], y=x[1], z=x[2]
dxdt[0] = -p*x[0] + p*x[1]; // dx/dt
dxdt[1] = -x[0] * x[2] + r*x[0] - x[1]; // dy/dt
dxdt[2] = x[0] * x[1] - b*x[2]; // dz/dt
return dxdt;
})
.set_integrate_range({ 0,25 })
.set_tolerance(8.0E-3L)
.set_step_size_range({ 1.0E-6,0.05 })
.integrate(state0,
[&](auto t, state_type state) {
timelog.emplace_back(t);
statelog.emplace_back(state);
});
return exit_status ? [&]{
// バイナリで出力
std::ofstream ofs(R"(./rkf45_test/lorenz.dat)", std::ios::binary);
const size_t N = timelog.size();
ofs.write(reinterpret_cast<const char*>(&N), sizeof(N));
ofs.write(reinterpret_cast<char*>(&timelog[0]), sizeof(long double)*N);
ofs.write(reinterpret_cast<char*>(&statelog[0]), sizeof(state_type)*N);
return EXIT_SUCCESS;
}()
: EXIT_FAILURE;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment