rk5adaptive.cc

Go to the documentation of this file.

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

#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")
  ;

  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);

  RK5ODEIntegrator<HarmonicOscillator> rk5(ho);

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

  cout << setprecision(16);

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

    cout << it.GetX() << ' '
         << it.GetY()[0] << ' '
         << it.GetY()[1] << '\n';

  return EXIT_SUCCESS;
}