//////////////////////////// // Anharmonic Oscillator // includes harmonic oscillator as a special case // // numerical integration the equations of motin using namespace std; #include // for cerr, cout #include // for setw() #include // for M_PI and cos() // acceleration = - (g/L) * x // also knowns as torque // assume units such that g/L equals 1 #define accel(x) (-(x)) const double pi(M_PI); void dumpit(double const period, int const Nst, int const steps_per_cycle, double const t, double const x, double const v, double const resid){ double denom(Nst); if (resid == 0. && denom == 0) denom = 1; // 0/0 is harmless; treat as 0/1 cout << fixed << setw(10) << t/period << scientific << setw(15) << x << scientific << setw(15) << v << scientific << setw(15) << (x*x + v*v -1.) << scientific << setw(15) << sqrt(resid / denom) << scientific << setw(15) << resid / denom << scientific << setw(15) << steps_per_cycle << endl; } int main(){ // principal configuration options: int verbosity(0); double Cmax(3); // desired number of cycles int steps_per_cycle(1000); // integrator steps per cycle // various quantities derived from the configuration: double period(2*pi); double dt(period/steps_per_cycle); // time per integrator step // initial conditions: double t(0); double x(1); // the position, also known as theta double v(0); // the velocity double resid(0); double Nst(0); // total number of integrator steps taken // prime the pump: double w(v + accel(x) * dt / 2.); // advanced velocity double oldw(w); // column headers: cout << setw(10) << "t/period" << setw(15) << "x " << setw(15) << "v " << setw(15) << " Delta E " << setw(15) << " rms resid " << setw(15) << " ms resid " << setw(15) << " steps/cycle" << endl; // outer loop over all cycles for (int cc = 0; cc < Cmax; cc++) { // inner loop; for explanation, see // https://www.av8n.com/physics/symplectic-integrator.htm#sec-verlet // // At the beg of the step, we have x[ii], t[ii], and w[ii] // At the end of the step, we have x[1+ii], t[1+ii], and w[1+ii] for (int ii = 0; ii < steps_per_cycle; ii++){ if (verbosity) { dumpit(period, Nst, steps_per_cycle, t, x, v, resid); } oldw = w; // leave some breadcrumbs in the forest x += w * dt; w += accel(x) * dt; t += dt; Nst++; double delta = x - cos(t); resid += delta*delta; // v is needed for printouts; otherwise wouldn't be needed: v = (oldw + w) / 2.; } // At this point, // oldw is w[Nsteps-1] // x is x[Nsteps] // w is w[Nsteps] v = (oldw + w) / 2.; if (cc == Cmax-1 || !verbosity) dumpit(period, Nst, steps_per_cycle, t, x, v, resid); } // end outer loop over cycles }