rkf45_test.cpp

#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;
}