Skip to content

Instantly share code, notes, and snippets.

@dc1394
Last active December 27, 2023 10:09
Show Gist options
  • Save dc1394/994e19418478fa14854b1a74d4f3df35 to your computer and use it in GitHub Desktop.
Save dc1394/994e19418478fa14854b1a74d4f3df35 to your computer and use it in GitHub Desktop.
Twitterの1D-Newton法の速度比較のコード(C++版、動的配列時間計測付き)
#include <chrono> // for std::chrono
#include <fstream> // for std::ofstream
#include <iomanip> // for std::setprecision
#include <ios> // for std::ios::fixed, std::ios::floatfield, std::scientific
#include <iostream> // for std::cout
#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()
{
using namespace std::chrono;
using clock = std::chrono::high_resolution_clock;
auto const mass = 1.0;
auto const k = 1.0;
auto const dt = 1e-2;
auto const nt = 100000000;
clock::time_point cp[5];
cp[0] = clock::now();
auto xt = new double[nt + 1];
auto vt = new double[nt + 1];
cp[1] = clock::now();
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;
}
cp[2] = clock::now();
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';
}
cp[3] = clock::now();
delete [] xt;
delete [] vt;
cp[4] = clock::now();
std::cout.setf(std::ios::fixed, std::ios::floatfield);
for (auto i = 0; i < 4; i++) {
std::cout << std::setprecision(4)
<< "Elapsed time cp[" << i << "] = "
<< duration_cast<duration<double>>(cp[i + 1] - cp[i]).count()
<< " (sec)\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;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment