ODEIntegrator.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <cmath>
#include <boost/program_options.hpp>

#include <utl/RK4ODEIntegrator.h>
#include <utl/RK5ODEIntegrator.h>

#include "HarmonicOscillator.h"

using namespace std;
namespace po = boost::program_options;
using namespace utl;


int
main(int argc, char* argv[])
{
  double dx;
  double accuracy;
  double x1;

  po::options_description desc("Allowed options");
  desc.add_options()
    ("help,h", "produce help message")
    ("dx,d", po::value<double>(&dx)->default_value(0.5), "step size")
    ("x1,1", po::value<double>(&x1)->default_value(1000), "end point")
    ("accuracy,a", po::value<double>(&accuracy)->default_value(1e-3),
     "variable step accuracy")
    ("variable,v", "run variable step")
    ("fixed,f", "run variable step with fixed try")
  ;

  po::variables_map vm;
  po::store(po::parse_command_line(argc, argv, desc), vm);
  po::notify(vm);

  if (vm.count("help")) {
    cout << desc << '\n';
    return EXIT_SUCCESS;
  }

  typedef double Vector[2];

  const double m = 1;
  const double k = 1;
  HarmonicOscillator ho(m, k);
  const double omega = sqrt(k / m);

  RK5ODEIntegrator<HarmonicOscillator> rk5(ho);

  const double x0 = 0;
  const Vector y0 = { 1, 0 };

  cout << setprecision(16);

  if (!vm.count("variable")) {

    RK4ODEIntegrator<HarmonicOscillator> rk4(ho);

    RK4Iterator<HarmonicOscillator, Vector> it4 = rk4.Begin(x0, y0);
    RK5Iterator<HarmonicOscillator, Vector> it5 = rk5.Begin(x0, y0);

    for ( ; it4.GetX() <= x1; it4 += dx, it5 += dx) {
      const double x = it4.GetX();
      cout << x << ' ' << it5.GetX() << ' '
           << it4.GetY()[0] << ' ' << it5.GetY()[0] << ' '
           << it4.GetY()[1] << ' ' << it5.GetY()[1] << ' '
           <<  y0[0] * cos(omega * x) << ' '
           << -y0[0] * sin(omega * x) << '\n';
    }

  } else {

    AdaptiveRK5Iterator<HarmonicOscillator, Vector> it =
      rk5.AdaptiveBegin(x0, dx, y0, accuracy);

    const bool fixed = vm.count("fixed");

    double x = it.GetX();

    while (x <= x1) {
      cout << x << ' '
           << it.GetY()[0] << ' '
           << it.GetY()[1] << ' '
           <<  y0[0] * cos(omega * x) << ' '
           << -y0[0] * sin(omega * x) << ' ';
      if (fixed)
        it += dx;
      else
        ++it;
      const double newx = it.GetX();
      cout << newx - x << '\n';
      x = newx;
    }

  }

  return EXIT_SUCCESS;
}