fork download
copy 
  1. #include <iostream>
  2. #include <cmath>
  3.  
  4. class equation {
  5. public:
  6.   double dxdt(double x, double y, double z, double vx, double vy, double vz) { return vx; }
  7.   double dydt(double x, double y, double z, double vx, double vy, double vz) { return vy; }
  8.   double dzdt(double x, double y, double z, double vx, double vy, double vz) { return vz; }
  9.   double dvxdt(double x, double y, double z, double vx, double vy, double vz) {
  10.     double r = std::sqrt(x * x + y * y + z * z);
  11.     return -x / (r * r * r);
  12.   }
  13.   double dvydt(double x, double y, double z, double vx, double vy, double vz) {
  14.     double r = std::sqrt(x * x + y * y + z * z);
  15.     return -y / (r * r * r);
  16.   }
  17.   double dvzdt(double x, double y, double z, double vx, double vy, double vz) {
  18.     double r = std::sqrt(x * x + y * y + z * z);
  19.     return -z / (r * r * r);
  20.   }
  21. };
  22.  
  23. struct mass {
  24.   double x, y, z, vx, vy, vz;
  25.   mass() {}
  26.   mass(double x, double y, double z, double vx, double vy, double vz) : x(x), y(y), z(z), vx(vx), vy(vy), vz(vz) {}
  27.   friend double area(mass m1, mass m2) {
  28.     double t1, t2, t3;
  29.     t1 = m1.y * m2.z - m1.z * m2.y;
  30.     t2 = m1.z * m2.x - m1.x * m2.z;
  31.     t3 = m1.x * m2.y - m1.y * m2.x;
  32.     return 0.5 * std::sqrt(t1 * t1 + t2 * t2 + t3 * t3);
  33.   }
  34.   friend std::ostream &operator<<(std::ostream &stream, mass m) {
  35.     stream << "(x, y, z, vx, vy, vz) = (" << m.x << ", " << m.y << ", " << m.z << ", ";
  36.     stream                                << m.vx << ", " << m.vy << ", " << m.vz << ")";
  37.     return stream;
  38.   }
  39. };
  40.  
  41. static int const N = 30000;
  42. static double const T  = 3.0;
  43.  
  44. class RungeKutta4 {
  45. public:
  46.   void operator()(mass &m, equation f) {
  47.     double k1x, k1y, k1z, k1vx, k1vy, k1vz;
  48.     double k2x, k2y, k2z, k2vx, k2vy, k2vz;
  49.     double k3x, k3y, k3z, k3vx, k3vy, k3vz;
  50.     double k4x, k4y, k4z, k4vx, k4vy, k4vz;
  51.     double h = T / N;
  52.  
  53.     k1x = h * f.dxdt(m.x, m.y, m.z, m.vx, m.vy, m.vz);
  54.     k1y = h * f.dydt(m.x, m.y, m.z, m.vx, m.vy, m.vz);
  55.     k1z = h * f.dzdt(m.x, m.y, m.z, m.vx, m.vy, m.vz);
  56.     k1vx = h * f.dvxdt(m.x, m.y, m.z, m.vx, m.vy, m.vz);
  57.     k1vy = h * f.dvydt(m.x, m.y, m.z, m.vx, m.vy, m.vz);
  58.     k1vz = h * f.dvzdt(m.x, m.y, m.z, m.vx, m.vy, m.vz);
  59.  
  60.     k2x = h * f.dxdt(m.x + k1x / 2.0, m.y + k1y / 2.0, m.z + k1z / 2.0,
  61.                      m.vx + k1vx / 2.0, m.vy + k1vy / 2.0, m.vz + k1vz / 2.0);
  62.     k2y = h * f.dydt(m.x + k1x / 2.0, m.y + k1y / 2.0, m.z + k1z / 2.0,
  63.                      m.vx + k1vx / 2.0, m.vy + k1vy / 2.0, m.vz + k1vz / 2.0);
  64.     k2z = h * f.dzdt(m.x + k1x / 2.0, m.y + k1y / 2.0, m.z + k1z / 2.0,
  65.                      m.vx + k1vx / 2.0, m.vy + k1vy / 2.0, m.vz + k1vz / 2.0);
  66.     k2vx = h * f.dvxdt(m.x + k1x / 2.0, m.y + k1y / 2.0, m.z + k1z / 2.0,
  67.                        m.vx + k1vx / 2.0, m.vy + k1vy / 2.0, m.vz + k1vz / 2.0);
  68.     k2vy = h * f.dvydt(m.x + k1x / 2.0, m.y + k1y / 2.0, m.z + k1z / 2.0,
  69.                        m.vx + k1vx / 2.0, m.vy + k1vy / 2.0, m.vz + k1vz / 2.0);
  70.     k2vz = h * f.dvzdt(m.x + k1x / 2.0, m.y + k1y / 2.0, m.z + k1z / 2.0,
  71.                        m.vx + k1vx / 2.0, m.vy + k1vy / 2.0, m.vz + k1vz / 2.0);
  72.  
  73.     k3x = h * f.dxdt(m.x + k2x / 2.0, m.y + k2y / 2.0, m.z + k2z,
  74.                      m.vx + k2vx / 2.0, m.vy + k2vy / 2.0, m.vz + k2vz / 2.0);
  75.     k3y = h * f.dydt(m.x + k2x / 2.0, m.y + k2y / 2.0, m.z + k2z,
  76.                      m.vx + k2vx / 2.0, m.vy + k2vy / 2.0, m.vz + k2vz / 2.0);
  77.     k3z = h * f.dzdt(m.x + k2x / 2.0, m.y + k2y / 2.0, m.z + k2z,
  78.                      m.vx + k2vx / 2.0, m.vy + k2vy / 2.0, m.vz + k2vz / 2.0);
  79.     k3vx = h * f.dvxdt(m.x + k2x / 2.0, m.y + k2y / 2.0, m.z + k2z,
  80.                        m.vx + k2vx / 2.0, m.vy + k2vy / 2.0, m.vz + k2vz / 2.0);
  81.     k3vy = h * f.dvydt(m.x + k2x / 2.0, m.y + k2y / 2.0, m.z + k2z,
  82.                        m.vx + k2vx / 2.0, m.vy + k2vy / 2.0, m.vz + k2vz / 2.0);
  83.     k3vz = h * f.dvzdt(m.x + k2x / 2.0, m.y + k2y / 2.0, m.z + k2z,
  84.                        m.vx + k2vx / 2.0, m.vy + k2vy / 2.0, m.vz + k2vz / 2.0);
  85.  
  86.     k4x = h * f.dxdt(m.x + k3x, m.y + k3y, m.z + k3z, m.vx + k3vx, m.vy + k3vy, m.vz + k3vz);
  87.     k4y = h * f.dydt(m.x + k3x, m.y + k3y, m.z + k3z, m.vx + k3vx, m.vy + k3vy, m.vz + k3vz);
  88.     k4z = h * f.dzdt(m.x + k3x, m.y + k3y, m.z + k3z, m.vx + k3vx, m.vy + k3vy, m.vz + k3vz);
  89.     k4vx = h * f.dvxdt(m.x + k3x, m.y + k3y, m.z + k3z, m.vx + k3vx, m.vy + k3vy, m.vz + k3vz);
  90.     k4vy = h * f.dvydt(m.x + k3x, m.y + k3y, m.z + k3z, m.vx + k3vx, m.vy + k3vy, m.vz + k3vz);
  91.     k4vz = h * f.dvzdt(m.x + k3x, m.y + k3y, m.z + k3z, m.vx + k3vx, m.vy + k3vy, m.vz + k3vz);
  92.  
  93.     m.x = m.x + (k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0;
  94.     m.y = m.y + (k1y + 2.0 * k2y + 2.0 * k3y + k4y) / 6.0;
  95.     m.z = m.z + (k1z + 2.0 * k2z + 2.0 * k3z + k4z) / 6.0;
  96.     m.vx = m.vx + (k1vx + 2.0 * k2vx + 2.0 * k3vx + k4vx) / 6.0;
  97.     m.vy = m.vy + (k1vy + 2.0 * k2vy + 2.0 * k3vy + k4vy) / 6.0;
  98.     m.vz = m.vz + (k1vz + 2.0 * k2vz + 2.0 * k3vz + k4vz) / 6.0;
  99.   }
  100. };
  101.  
  102. int const K = 2000;
  103. int main() {
  104.   RungeKutta4 e;
  105.  
  106.   equation f;
  107.   mass m(0.0, 0.70710678, 0.70710678, 1.0, 0.0, 0.0), m0;
  108.   std::cout << "t = 0.0 : " << m << std::endl;
  109.   int i = K + 1;
  110.   double s = 0.0;
  111.   for (int k = 0; k < T * N; k++) {
  112.     m0 = m;
  113.     e(/*ref*/m, f);
  114.     s += area(m0, m);
  115.     if (--i == 0) {
  116.       std::cout << "t = " << (double)k / N << ':' << m << ": s = " << s << std::endl;
  117.       i = K;
  118.       s = 0.0;
  119.     }
  120.   }
  121.   return 0;
  122. }
  123. /* end */
  124.  
Success #stdin #stdout 0.03s 2856KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
t = 0.0 : (x, y, z, vx, vy, vz) = (0, 0.707107, 0.707107, 1, 0, 0)
t = 0.0666667:(x, y, z, vx, vy, vz) = (0.198767, 0.692998, 0.692998, 0.980047, -0.14055, -0.14055): s = 0.10005
t = 0.133333:(x, y, z, vx, vy, vz) = (0.38951, 0.651261, 0.651261, 0.921022, -0.275427, -0.275426): s = 0.1
t = 0.2:(x, y, z, vx, vy, vz) = (0.564725, 0.58356, 0.58356, 0.825278, -0.399324, -0.399322): s = 0.1
t = 0.266667:(x, y, z, vx, vy, vz) = (0.717425, 0.492594, 0.492595, 0.696632, -0.507301, -0.507297): s = 0.0999999
t = 0.333333:(x, y, z, vx, vy, vz) = (0.841523, 0.38199, 0.381992, 0.540212, -0.595053, -0.595048): s = 0.0999998
t = 0.4:(x, y, z, vx, vy, vz) = (0.932072, 0.256156, 0.256159, 0.362256, -0.659083, -0.659076): s = 0.0999997
t = 0.466667:(x, y, z, vx, vy, vz) = (0.985461, 0.120111, 0.120115, 0.169857, -0.696836, -0.696828): s = 0.0999996
t = 0.533333:(x, y, z, vx, vy, vz) = (0.999563, -0.0207235, -0.0207174, -0.0293147, -0.706808, -0.706799): s = 0.0999994
t = 0.6:(x, y, z, vx, vy, vz) = (0.973814, -0.160732, -0.160723, -0.227317, -0.688602, -0.688591): s = 0.0999992
t = 0.666667:(x, y, z, vx, vy, vz) = (0.909241, -0.294332, -0.294321, -0.416258, -0.642942, -0.642931): s = 0.0999991
t = 0.733333:(x, y, z, vx, vy, vz) = (0.808418, -0.416197, -0.416185, -0.588604, -0.571649, -0.571638): s = 0.0999989
t = 0.8:(x, y, z, vx, vy, vz) = (0.675365, -0.52147, -0.521455, -0.737485, -0.477565, -0.477555): s = 0.0999988
t = 0.866667:(x, y, z, vx, vy, vz) = (0.515387, -0.605952, -0.605936, -0.856964, -0.36444, -0.364433): s = 0.0999988
t = 0.933333:(x, y, z, vx, vy, vz) = (0.334861, -0.666276, -0.666258, -0.942278, -0.236785, -0.236779): s = 0.0999987
t = 1:(x, y, z, vx, vy, vz) = (0.140984, -0.700036, -0.700017, -0.990024, -0.0996876, -0.0996851): s = 0.0999987
t = 1.06667:(x, y, z, vx, vy, vz) = (-0.0585141, -0.705885, -0.705867, -0.9983, 0.0413858, 0.0413847): s = 0.0999987
t = 1.13333:(x, y, z, vx, vy, vz) = (-0.255679, -0.683591, -0.683574, -0.966772, 0.180811, 0.180806): s = 0.0999987
t = 1.2:(x, y, z, vx, vy, vz) = (-0.44265, -0.634042, -0.634026, -0.896699, 0.313028, 0.313019): s = 0.0999987
t = 1.26667:(x, y, z, vx, vy, vz) = (-0.611972, -0.559213, -0.559199, -0.790872, 0.432766, 0.432752): s = 0.0999986
t = 1.33333:(x, y, z, vx, vy, vz) = (-0.756894, -0.462088, -0.462077, -0.653511, 0.53525, 0.535232): s = 0.0999986
t = 1.4:(x, y, z, vx, vy, vz) = (-0.871637, -0.346538, -0.346531, -0.490093, 0.616393, 0.616372): s = 0.0999985
t = 1.46667:(x, y, z, vx, vy, vz) = (-0.951626, -0.217171, -0.217169, -0.307133, 0.67296, 0.672936): s = 0.0999984
t = 1.53333:(x, y, z, vx, vy, vz) = (-0.993672, -0.0791449, -0.0791477, -0.111926, 0.702695, 0.702668): s = 0.0999982
t = 1.6:(x, y, z, vx, vy, vz) = (-0.996099, 0.0620367, 0.0620284, 0.0877441, 0.704412, 0.704383): s = 0.099998
t = 1.66667:(x, y, z, vx, vy, vz) = (-0.958808, 0.200745, 0.200731, 0.283916, 0.678042, 0.678014): s = 0.0999979
t = 1.73333:(x, y, z, vx, vy, vz) = (-0.883288, 0.331449, 0.331429, 0.468767, 0.624637, 0.62461): s = 0.0999977
t = 1.8:(x, y, z, vx, vy, vz) = (-0.772549, 0.448937, 0.448912, 0.634928, 0.546326, 0.546302): s = 0.0999976
t = 1.86667:(x, y, z, vx, vy, vz) = (-0.631007, 0.548525, 0.548496, 0.775773, 0.446232, 0.446211): s = 0.0999975
t = 1.93333:(x, y, z, vx, vy, vz) = (-0.464304, 0.626242, 0.626209, 0.885687, 0.328345, 0.328329): s = 0.0999974
t = 2:(x, y, z, vx, vy, vz) = (-0.279089, 0.678989, 0.678954, 0.960286, 0.197366, 0.197356): s = 0.0999974
t = 2.06667:(x, y, z, vx, vy, vz) = (-0.082746, 0.704663, 0.704626, 0.996596, 0.0585163, 0.0585133): s = 0.0999974
t = 2.13333:(x, y, z, vx, vy, vz) = (0.116897, 0.70224, 0.702203, 0.993169, -0.0826672, -0.0826628): s = 0.0999974
t = 2.2:(x, y, z, vx, vy, vz) = (0.311878, 0.671816, 0.671781, 0.950141, -0.220555, -0.220543): s = 0.0999974
t = 2.26667:(x, y, z, vx, vy, vz) = (0.494424, 0.614605, 0.614573, 0.869228, -0.349649, -0.34963): s = 0.0999974
t = 2.33333:(x, y, z, vx, vy, vz) = (0.657256, 0.532887, 0.53286, 0.753655, -0.464802, -0.464775): s = 0.0999973
t = 2.4:(x, y, z, vx, vy, vz) = (0.79388, 0.429921, 0.4299, 0.60803, -0.561422, -0.561389): s = 0.0999972
t = 2.46667:(x, y, z, vx, vy, vz) = (0.898849, 0.309813, 0.309799, 0.438161, -0.635656, -0.635617): s = 0.0999971
t = 2.53333:(x, y, z, vx, vy, vz) = (0.967978, 0.177351, 0.177345, 0.25082, -0.684544, -0.684501): s = 0.099997
t = 2.6:(x, y, z, vx, vy, vz) = (0.998509, 0.0378171, 0.0378205, 0.0534776, -0.706137, -0.706091): s = 0.0999969
t = 2.66667:(x, y, z, vx, vy, vz) = (0.989226, -0.103225, -0.103212, -0.145997, -0.699573, -0.699527): s = 0.0999967
t = 2.73333:(x, y, z, vx, vy, vz) = (0.940499, -0.24015, -0.240128, -0.339651, -0.665115, -0.665069): s = 0.0999965
t = 2.8:(x, y, z, vx, vy, vz) = (0.85427, -0.3675, -0.367469, -0.519762, -0.604135, -0.604092): s = 0.0999964
t = 2.86667:(x, y, z, vx, vy, vz) = (0.733978, -0.480196, -0.480157, -0.679149, -0.519065, -0.519028): s = 0.0999963
t = 2.93333:(x, y, z, vx, vy, vz) = (0.584418, -0.573744, -0.573699, -0.811456, -0.413297, -0.413267): s = 0.0999962
https://ideone.com/E6LEsO
language:
C++ (gcc 8.3)
created:
11 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!