Twitterの1D-Newton法の速度比較のコード

newton.cpp

#include <fstream>  // for std::ofstream
#include <iomanip>  // for std::setprecision
#include <ios>      // for std::scientific
#include <utility>  // for std::make_pair, std::pair
#include <vector>   // for std::vector

namespace {
    std::pair<double, double> Runge_Kutta_4th(double x, double v, double dt, double mass, double k);
    double force(double x, double mass, double k);
}

int main()
{
    auto const mass = 1.0;
    auto const k = 1.0;
    auto const dt = 1e-2;
    auto const nt = 100000000;

    std::vector<double> xt(nt + 1);
    std::vector<double> vt(nt + 1);

    auto x = 0.0;
    auto v = 1.0;

    for (auto it = 0; it <= nt; it++) {
        xt[it] = x;
        vt[it] = v;
        auto const ret = Runge_Kutta_4th(x, v, dt, mass, k);
        x = ret.first;
        v = ret.second;
    }

    std::ofstream ofs("result_cpp.out");
    ofs.setf(std::ios::scientific);
    for (auto it = nt - 1000; it <= nt; it++) {
        ofs << std::setprecision(7) << static_cast<double>(it) * dt << ' ' << xt[it] << ' ' << vt[it] << '\n';
    }

    return 0;
}

namespace {
    std::pair<double, double> Runge_Kutta_4th(double x, double v, double dt, double mass, double k)
    {
        auto const x1 = v;
        auto const v1 = force(x, mass, k);

        auto const x2 = v + 0.5 * dt * v1;
        auto const v2 = force(x + 0.5 * x1 * dt, mass, k);

        auto const x3 = v + 0.5 * dt * v2;
        auto const v3 = force(x + 0.5 * x2 * dt, mass, k);

        auto const x4 = v + dt * v3;
        auto const v4 = force(x + x3 * dt, mass, k);

        x += (x1 + 2 * x2 + 2 * x3 + x4) * dt / 6.0;
        v += (v1 + 2 * v2 + 2 * v3 + v4) * dt / 6.0;

        return std::make_pair(x, v);
    }

    double force(double x, double mass, double k)
    {
        return -x * k / mass;
    }
}