fork download 
  1. #include <iostream>
  2. #include <cmath>
  3.  
  4. class equation {
  5. public:
  6.   double dxdt(double x, double y, double vx, double vy) { return vx; }
  7.   double dydt(double x, double y, double vx, double vy) { return vy; }
  8.   double dvxdt(double x, double y, double vx, double vy) {
  9.     double r = sqrt(x * x + y * y);
  10.     return -x / (r * r * r);
  11.   }
  12.   double dvydt(double x, double y, double u, double v) {
  13.     double r = sqrt(x * x + y * y);
  14.     return -y / (r * r * r);
  15.   }
  16. };
  17.  
  18. struct mass {
  19.   double x, y, vx, vy;
  20.   mass() {}
  21.   mass(double x, double y, double vx, double vy) : x(x), y(y), vx(vx), vy(vy) {}
  22.   friend double area(mass m1, mass m2) {
  23.     return 0.5 * std::fabs(m1.x * m2.y - m2.x * m1.y);
  24.   }
  25.   friend std::ostream &operator<<(std::ostream &stream, mass m) {
  26.     stream.width(9);
  27.     stream.precision(6);
  28.     stream << "(x, y, vx, vy) = (" << m.x << ", " << m.y << ", ";
  29.     stream                         << m.vx << ", " << m.vy << ")";
  30.     return stream;
  31.   }
  32. };
  33.  
  34. static int const N = 10000;
  35. static double const T  = 10.0;
  36.  
  37. class RungeKutta2 {
  38.   double k1x, k2x, k1y, k2y, k1vx, k2vx, k1vy, k2vy;
  39. public:
  40.   void operator()(mass &m, equation f) {
  41.     double h = T / N;
  42.     k1x = f.dxdt(m.x, m.y, m.vx, m.vy);
  43.     k1y = f.dydt(m.x, m.y, m.vx, m.vy);
  44.     k1vx = f.dvxdt(m.x, m.y, m.vx, m.vy);
  45.     k1vy = f.dvydt(m.x, m.y, m.vx, m.vy);
  46.  
  47.     k2x = f.dxdt(m.x + h * k1x, m.y + h * k1y, m.vx + h * k1vx, m.vy + h * k1vy);
  48.     k2y = f.dydt(m.x + h * k1x, m.y + h * k1y, m.vx + h * k1vy, m.vy + h * k1vy);
  49.     k2vx = f.dvxdt(m.x + h * k1x, m.y + h * k1y, m.vx + h * k1vx, m.vy + h * k1vy);
  50.     k2vy = f.dvydt(m.x + h * k1x, m.y + h * k1y, m.vx + h * k1vx, m.vy + h * k1vy);
  51.  
  52.     m.x = m.x + h * (k1x + k2x) / 2.0;
  53.     m.y = m.y + h * (k1y + k2y) / 2.0;
  54.     m.vx = m.vx + h * (k1vx + k2vx) / 2.0;
  55.     m.vy = m.vy + h * (k1vy + k2vy) / 2.0;
  56.   }
  57. };
  58.  
  59. int const K = 1000;
  60. int main() {
  61.   RungeKutta2 euler;
  62.  
  63.   equation f;
  64.   mass m(0.0, 1.0, 1.2, 0.0), m0;
  65.   int i = K + 1;
  66.   double s = 0.0;
  67.   for (int k = 0; k < T * N; k++) {
  68.     m0 = m;
  69.     euler(/*ref*/m, f);
  70.     s += area(m0, m);
  71.     if (--i == 0) {
  72.       std::cout << "t = " << (double)k / N << ':' << m << ": s = " << s << std::endl;
  73.       i = K;
  74.       s = 0.0;
  75.     }
  76.   }
  77.   return 0;
  78. }
  79. /* end */
  80.  
Success #stdin #stdout 0.01s 2852KB
comments (?)
stdin 
Standard input is empty
stdout 
t = 0.1:(x, y, vx, vy) = (1.03847, 0.574968, 0.770319, -0.729047): s = 0.6006
t = 0.2:(x, y, vx, vy) = (1.5253, -0.234641, 0.239964, -0.823644): s = 0.6
t = 0.3:(x, y, vx, vy) = (1.59136, -1.00564, -0.0785058, -0.704461): s = 0.6
t = 0.4:(x, y, vx, vy) = (1.41231, -1.63148, -0.263387, -0.545413): s = 0.6
t = 0.5:(x, y, vx, vy) = (1.08944, -2.09605, -0.372753, -0.384319): s = 0.6
t = 0.6:(x, y, vx, vy) = (0.682371, -2.40169, -0.434939, -0.227753): s = 0.6
t = 0.7:(x, y, vx, vy) = (0.230756, -2.55285, -0.463282, -0.0750202): s = 0.6
t = 0.8:(x, y, vx, vy) = (-0.234764, -2.55219, -0.463163, 0.0763324): s = 0.6
t = 0.9:(x, y, vx, vy) = (-0.686132, -2.39971, -0.434559, 0.229089): s = 0.6
t = 1:(x, y, vx, vy) = (-1.09266, -2.09271, -0.372037, 0.385696): s = 0.6
t = 1.1:(x, y, vx, vy) = (-1.41458, -1.62675, -0.262168, 0.546819): s = 0.6
t = 1.2:(x, y, vx, vy) = (-1.59203, -0.999531, -0.0764358, 0.705765): s = 0.6
t = 1.3:(x, y, vx, vy) = (-1.52321, -0.227506, 0.243566, 0.824191): s = 0.6
t = 1.4:(x, y, vx, vy) = (-1.03177, 0.581269, 0.775698, 0.726043): s = 0.6
t = 1.5:(x, y, vx, vy) = (0.00919407, 0.999971, 1.19996, -0.00766081): s = 0.6
t = 1.6:(x, y, vx, vy) = (1.04359, 0.570105, 0.766184, -0.73132): s = 0.6
t = 1.7:(x, y, vx, vy) = (1.52689, -0.240124, 0.237205, -0.823214): s = 0.6
t = 1.8:(x, y, vx, vy) = (1.59084, -1.01032, -0.080091, -0.703456): s = 0.6
t = 1.9:(x, y, vx, vy) = (1.41055, -1.6351, -0.264321, -0.544332): s = 0.6
t = 2:(x, y, vx, vy) = (1.08696, -2.0986, -0.373302, -0.383261): s = 0.6
t = 2.1:(x, y, vx, vy) = (0.679477, -2.4032, -0.43523, -0.226726): s = 0.6
t = 2.2:(x, y, vx, vy) = (0.227674, -2.55334, -0.463373, -0.0740111): s = 0.6
t = 2.3:(x, y, vx, vy) = (-0.237844, -2.55168, -0.463069, 0.0773416): s = 0.6
t = 2.4:(x, y, vx, vy) = (-0.689021, -2.39818, -0.434264, 0.230117): s = 0.6
t = 2.5:(x, y, vx, vy) = (-1.09513, -2.09014, -0.371483, 0.386755): s = 0.6
t = 2.6:(x, y, vx, vy) = (-1.41632, -1.6231, -0.261226, 0.547901): s = 0.6
t = 2.7:(x, y, vx, vy) = (-1.59253, -0.994828, -0.0748374, 0.706767): s = 0.6
t = 2.8:(x, y, vx, vy) = (-1.52157, -0.222016, 0.246348, 0.824603): s = 0.6
t = 2.9:(x, y, vx, vy) = (-1.02659, 0.586097, 0.779837, 0.723697): s = 0.6
t = 3:(x, y, vx, vy) = (0.0171877, 0.999897, 1.19988, -0.0143208): s = 0.6
t = 3.1:(x, y, vx, vy) = (1.04868, 0.565227, 0.762051, -0.733562): s = 0.6
t = 3.2:(x, y, vx, vy) = (1.52846, -0.245605, 0.234457, -0.822777): s = 0.6
t = 3.3:(x, y, vx, vy) = (1.5903, -1.015, -0.0816704, -0.702449): s = 0.6
t = 3.4:(x, y, vx, vy) = (1.40879, -1.63872, -0.265251, -0.54325): s = 0.6
t = 3.5:(x, y, vx, vy) = (1.08447, -2.10115, -0.373849, -0.382202): s = 0.6
t = 3.6:(x, y, vx, vy) = (0.676581, -2.40471, -0.43552, -0.225699): s = 0.6
t = 3.7:(x, y, vx, vy) = (0.224592, -2.55383, -0.463462, -0.0730022): s = 0.6
t = 3.8:(x, y, vx, vy) = (-0.240924, -2.55116, -0.462975, 0.0783509): s = 0.6
t = 3.9:(x, y, vx, vy) = (-0.691909, -2.39665, -0.433969, 0.231144): s = 0.6
t = 4:(x, y, vx, vy) = (-1.0976, -2.08756, -0.370927, 0.387815): s = 0.6
t = 4.1:(x, y, vx, vy) = (-1.41805, -1.61945, -0.260282, 0.548982): s = 0.6
t = 4.2:(x, y, vx, vy) = (-1.59302, -0.990118, -0.0732332, 0.707768): s = 0.6
t = 4.3:(x, y, vx, vy) = (-1.51992, -0.216523, 0.24914, 0.825006): s = 0.6
t = 4.4:(x, y, vx, vy) = (-1.02138, 0.59091, 0.783977, 0.721318): s = 0.6
t = 4.5:(x, y, vx, vy) = (0.0251804, 0.99978, 1.19974, -0.0209792): s = 0.6
t = 4.6:(x, y, vx, vy) = (1.05374, 0.560334, 0.757921, -0.735772): s = 0.6
t = 4.7:(x, y, vx, vy) = (1.53002, -0.251083, 0.231718, -0.822331): s = 0.6
t = 4.8:(x, y, vx, vy) = (1.58975, -1.01968, -0.0832441, -0.701442): s = 0.6
t = 4.9:(x, y, vx, vy) = (1.40702, -1.64234, -0.266179, -0.542169): s = 0.6
t = 5:(x, y, vx, vy) = (1.08198, -2.10369, -0.374394, -0.381144): s = 0.6
t = 5.1:(x, y, vx, vy) = (0.673683, -2.40621, -0.435808, -0.224672): s = 0.6
t = 5.2:(x, y, vx, vy) = (0.221509, -2.55431, -0.46355, -0.0719933): s = 0.6
t = 5.3:(x, y, vx, vy) = (-0.244003, -2.55064, -0.462879, 0.0793602): s = 0.6
t = 5.4:(x, y, vx, vy) = (-0.694794, -2.3951, -0.433671, 0.232172): s = 0.6
t = 5.5:(x, y, vx, vy) = (-1.10007, -2.08498, -0.37037, 0.388874): s = 0.6
t = 5.6:(x, y, vx, vy) = (-1.41978, -1.61579, -0.259334, 0.550064): s = 0.6
t = 5.7:(x, y, vx, vy) = (-1.5935, -0.985402, -0.0716231, 0.708767): s = 0.6
t = 5.8:(x, y, vx, vy) = (-1.51825, -0.211028, 0.251942, 0.825401): s = 0.6
t = 5.9:(x, y, vx, vy) = (-1.01614, 0.595706, 0.788118, 0.718907): s = 0.6
t = 6:(x, y, vx, vy) = (0.033172, 0.999618, 1.19954, -0.0276354): s = 0.6
t = 6.1:(x, y, vx, vy) = (1.05877, 0.555426, 0.753793, -0.737951): s = 0.6
t = 6.2:(x, y, vx, vy) = (1.53155, -0.256558, 0.22899, -0.821878): s = 0.6
t = 6.3:(x, y, vx, vy) = (1.58919, -1.02435, -0.0848122, -0.700432): s = 0.6
t = 6.4:(x, y, vx, vy) = (1.40525, -1.64594, -0.267103, -0.541088): s = 0.6
t = 6.5:(x, y, vx, vy) = (1.07949, -2.10622, -0.374937, -0.380086): s = 0.6
t = 6.6:(x, y, vx, vy) = (0.670783, -2.4077, -0.436094, -0.223645): s = 0.6
t = 6.7:(x, y, vx, vy) = (0.218425, -2.55479, -0.463637, -0.0709844): s = 0.6
t = 6.8:(x, y, vx, vy) = (-0.247081, -2.5501, -0.462782, 0.0803695): s = 0.6
t = 6.9:(x, y, vx, vy) = (-0.697678, -2.39356, -0.433373, 0.233201): s = 0.6
t = 7:(x, y, vx, vy) = (-1.10253, -2.08239, -0.36981, 0.389934): s = 0.6
t = 7.1:(x, y, vx, vy) = (-1.4215, -1.61213, -0.258383, 0.551145): s = 0.6
t = 7.2:(x, y, vx, vy) = (-1.59397, -0.98068, -0.0700073, 0.709764): s = 0.6
t = 7.3:(x, y, vx, vy) = (-1.51656, -0.20553, 0.254754, 0.825788): s = 0.6
t = 7.4:(x, y, vx, vy) = (-1.01088, 0.600486, 0.79226, 0.716463): s = 0.6
t = 7.5:(x, y, vx, vy) = (0.0411621, 0.999412, 1.19929, -0.0342888): s = 0.6
t = 7.6:(x, y, vx, vy) = (1.06378, 0.550504, 0.749669, -0.7401): s = 0.6
t = 7.7:(x, y, vx, vy) = (1.53307, -0.262029, 0.226271, -0.821417): s = 0.6
t = 7.8:(x, y, vx, vy) = (1.58863, -1.02901, -0.0863745, -0.699422): s = 0.6
t = 7.9:(x, y, vx, vy) = (1.40347, -1.64954, -0.268024, -0.540006): s = 0.6
t = 8:(x, y, vx, vy) = (1.077, -2.10875, -0.375477, -0.379028): s = 0.6
t = 8.1:(x, y, vx, vy) = (0.667882, -2.40918, -0.436379, -0.222619): s = 0.6
t = 8.2:(x, y, vx, vy) = (0.215341, -2.55526, -0.463723, -0.0699757): s = 0.6
t = 8.3:(x, y, vx, vy) = (-0.250159, -2.54957, -0.462684, 0.0813789): s = 0.6
t = 8.4:(x, y, vx, vy) = (-0.700559, -2.392, -0.433072, 0.234229): s = 0.6
t = 8.5:(x, y, vx, vy) = (-1.10498, -2.07979, -0.369248, 0.390994): s = 0.6
t = 8.6:(x, y, vx, vy) = (-1.42322, -1.60845, -0.257429, 0.552226): s = 0.6
t = 8.7:(x, y, vx, vy) = (-1.59443, -0.975951, -0.0683855, 0.71076): s = 0.6
t = 8.8:(x, y, vx, vy) = (-1.51486, -0.200029, 0.257576, 0.826166): s = 0.6
t = 8.9:(x, y, vx, vy) = (-1.00559, 0.605249, 0.796403, 0.713987): s = 0.6
t = 9:(x, y, vx, vy) = (0.0491502, 0.999161, 1.19899, -0.0409385): s = 0.6
t = 9.1:(x, y, vx, vy) = (1.06876, 0.545567, 0.745547, -0.742217): s = 0.6
t = 9.2:(x, y, vx, vy) = (1.53457, -0.267498, 0.223563, -0.820949): s = 0.6
t = 9.3:(x, y, vx, vy) = (1.58805, -1.03366, -0.0879312, -0.698411): s = 0.6
t = 9.4:(x, y, vx, vy) = (1.40169, -1.65313, -0.268942, -0.538925): s = 0.6
t = 9.5:(x, y, vx, vy) = (1.0745, -2.11127, -0.376016, -0.377971): s = 0.6
t = 9.6:(x, y, vx, vy) = (0.664978, -2.41066, -0.436663, -0.221592): s = 0.6
t = 9.7:(x, y, vx, vy) = (0.212257, -2.55572, -0.463807, -0.068967): s = 0.6
t = 9.8:(x, y, vx, vy) = (-0.253236, -2.54902, -0.462584, 0.0823884): s = 0.6
t = 9.9:(x, y, vx, vy) = (-0.703438, -2.39044, -0.432771, 0.235257): s = 0.6
https://ideone.com/jLJydK
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!