gistfile1.txt

#include <iostream>
#include <iomanip>
#include <valarray>
#include <iostream>
#include <fstream>
using namespace std;

// Function prototypes
void step( double dx, double &x, valarray<double> &Y );
valarray<double> F( double x, valarray<double> Y );

//======================================================================


int main(){
   int nstep;
   double x, dx;

   // Set number of variables
   int ndep = 4;
   valarray<double> Y(ndep);

   // Set initial values
   x = 0;   
   Y[0] = 1.0;
   Y[1] = 1.0;
   Y[2] = 0.1;
   Y[3] = 0;

   // Set step size and number of steps
   dx = 0.01;
   nstep = 1000;


   // Write header and initial conditions
   #define SP << setw( 12 ) << fixed << setprecision( 4 ) <<    // save some repetition when writing
   cout << 'x' SP 'u'  SP 'v'  << endl;
   cout <<  x  SP Y[0] SP Y[1] SP Y[2] SP Y[3]<< endl;

   // Solve the differential equation using nstep intervals
   ofstream myfile;
   myfile.open ("data.js");
   myfile << "let data = [\n";
   for ( int n = 1; n <= nstep; n++ ) {
      step( dx, x, Y );                          // main Runge-Kutta step
      cout << x SP Y[0] SP Y[1] SP Y[2] SP Y[3]<< endl;         // output
      myfile << "{\"x\": \"" << Y[0] << "\",\n\"y\": \"" << Y[1] << "\",\n\"vx\": \"" << Y[2] << "\",\n\"vy\": \"" << Y[3] << "\"\n},";
   }
   myfile.close();
}

//======================================================================

void step( double dx, double &x, valarray<double> &Y )
{
   int ndep = Y.size();
   valarray<double> dY1(ndep), dY2(ndep), dY3(ndep), dY4(ndep);

   dY1 = F( x           , Y             ) * dx;
   dY2 = F( x + 0.5 * dx, Y + 0.5 * dY1 ) * dx;
   dY3 = F( x + 0.5 * dx, Y + 0.5 * dY2 ) * dx;
   dY4 = F( x +       dx, Y       + dY3 ) * dx;
   Y += ( dY1 + 2.0 * dY2 + 2.0 * dY3 + dY4 ) / 6.0;

   x += dx;
}

//======================================================================

valarray<double> F( double x, valarray<double> Y ){
   valarray<double> f( Y.size() );         

   f[0] = Y[2];
   f[1] = Y[3];
   f[2] = -Y[0];
   f[3] = -Y[1];   
                                                
   return f;
}