desystest.cpp

Go to the documentation of this file.

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

#include <commandline.h>
#include <desys.h>
#include <desystest.h>
#include <desystestspring2.h>
#include <nintegration.h>
#include <typedefs.h>

// <TODO> Macro solution for 1D de's.
// ie design a macro to plug in equation and get class.

string desystest::doc[] =
{
  "",
  "Euler on y''+y=0 where y(0)=1 and y'(0)=0.", 
  "TrapezoidEuler on y''+y=0 where y(0)=1 and y'(0)=0.",
  "Point in time - RungKutta on y''+y=0 where y(0)=1 and y'(0)=0.", 
  "Point in time - Euler on y''+y=0 as 2 equations using desys. ",
  "Spring mass system with 2 masses. " 
};


template< typename T >
class desystestf3
{
public:

  T * xval;
  T * yi;

  desystestf3(T const * xval_, T const * yi_)
    : xval(xval_), yi(yi_) { *xval=0; yi[0]=1; }

  void operator ++ () {}

  void operator () (T &  Dz, T const z, T const x) const
    { Dz = z; }
};




void desystestprobleminfo()
{
  cout << "y''+y = 0" << endl;
  cout << "y(0)=1, y'(0)=0" << endl;
  cout << "y=cos(x) by algebra." << endl;
  cout << endl;
}

/*
template< class D >
void desystest00( D & de )
{
  de.xval = 0.0;

  de.yiGet(0,0) = 1.0;
  de.yiGet(1,0) = 0.0;

  desystestprobleminfo();

  cout << "Enter h: ";
  cin >> de.hi[0];

  double xmax;
  cout << "Enter xmax: ";
  cin >> xmax;

  for ( ; de.xval<xmax; )
    de();

  cout << SHOW(de.xval) << endl;
  cout << SHOW(de.yiGet(0,0)) << endl;
  cout << SHOW(de.yiGet(1,0)) << endl;
}
*/

void desystest::test01(int argc, char** argv)
{
  commandline cmd(argc,argv); 

  //typedef nintegrationRK< desystestf1<double>, double > integrator;
  typedef nintegrationRK< f1<double>, double > integrator;
  //typedef nintegrationRK< desystestf3<double>, double > integrator;
  
  //desys< eul ,double > de(1,2);
  desys< integrator ,double > de(1,2);

  //de.yiGet(0,0)=1.0;
  //de.yiGet(0,1)=0.0;

//  de.yiGet(0,0)=1.0;
//  de.yiGet(0,1)=0.0;
  de.hstep=0.1;
  cmd.mapvar(de.hstep,"h");
  double xmax(0.3);
  cmd.mapvar(xmax,"xmax");

  cout << SHOW(de.xval) << endl;
  cout << SHOW(de.hstep) << endl;
  cout << SHOW(xmax) << endl;

  //cout << de << endl;

  double delta=1e-10;
  for (; de.xval+delta < xmax; de.eval() )
  {
cout << "***---------***" << endl;
    cout << de << endl;
  }
  cout << de << endl;
}

void desystest::test02()
{
  typedef nintegrationTrapezoidEuler< f1<double>, double > eul;
  
  desys< eul ,double > de(1,2);

  de.yiGet(0,0)=1.0;
  de.yiGet(0,1)=0.0;
  de.hstep = 0.001;

  cout << de << endl;

  for (; de.xval < 0.3; de.eval() );
    cout << de << endl;
  cout << de << endl;
}

int desystest::test03(int argc, char** argv)
{
  typedef nintegrationRK< f1<double>, double > integrator;
  
  desys< integrator ,double > de(1,2);

  //de.yiGet(0,0)=1.0;
  //de.yiGet(0,1)=0.0;
  de.hstep=0.001;

  commandline cmd(argc,argv); 
  cmd.mapvar(de.hstep,"hstep");


  cout << de << endl;

  double time_end=0.3;
  cmd.mapvar(time_end,"time_end");
  cout << SHOW(time_end) << endl;

  // Iterate in time to time_end;
  for (; de.xval < time_end; de.eval() );

  cout << de << endl;

  double error[2];
  error[0] = abs(cos(de.xval) - de.yiGet(0,0));
  error[1] = abs(-sin(de.xval) - de.yiGet(0,1));

  cout << SHOW(error[0]) << endl;
  cout << SHOW(error[1]) << endl;

  double myerror=de.hstep*10;
  assertreturnOS( error[0] < myerror );
  assertreturnOS( error[1] < myerror );

  return 0;
}


class desystestf2
{
public:

  doublec * xval;
  doublec * yi;

  void equ01 (double & Dminus1, doublec y, doublec x ) const
    { Dminus1 = -yi[1]; }
  void equ02 (double & Dminus1, doublec y, doublec x ) const
    { Dminus1 = yi[0]; }

  uint current;

  desystestf2(doublec * xval_, doublec * yi_)
    : xval(xval_), yi(yi_), current(0)
    {}

  void operator () (double & Dminus1, doublec y, doublec x ) const
    { 
      switch (current)
      {
        case 0: equ01(Dminus1,y,x); break;
        case 1: equ02(Dminus1,y,x); break;
      }
    }

  void operator ++ ()
    { ++current; current %= 2; }
};


int desystest::test04(int argc, char** argv)
{

  typedef nintegrationEuler< desystestf2, double > eul;
  
  desys< eul ,double > de(2,1);

  de.yiGet(0,0)=0.0;
  de.yiGet(1,0)=1.0;
  de.hstep=0.00001;

  commandline cmd(argc,argv); 

  cmd.mapvar(de.hstep,"hstep");

  double time_end=0.3;
  cmd.mapvar(time_end,"time_end");
  cout << SHOW(time_end) << endl;

  cout << de << endl;

  for (; de.xval < time_end; de.eval() );

  cout << de << endl;


  double error[2];
  error[0] = abs(cos(de.xval) - de.yiGet(1,0));
  error[1] = abs(-sin(de.xval) - de.yiGet(0,0));

  cout << SHOW(error[0]) << endl;
  cout << SHOW(error[1]) << endl;

  double myerror=de.hstep*100;
  assertreturnOS( error[0] < myerror );
  assertreturnOS( error[1] < myerror );

  return 0;
}


void desystest::test05(int argc, char** argv)
{
  typedef nintegrationRK< desystestspring2, double > integrator;
  desys< integrator, double > de(2,2);

  commandline cmd(argc,argv); 

  double w0=1.0;
  double w1=0.0;
  double Dw0=0.3;
  double Dw1=0.2;
  double h=0.001;

  cmd.mapvar(w0,"w0");
  cmd.mapvar(Dw0,"Dw0");
  cmd.mapvar(w1,"w1");
  cmd.mapvar(Dw1,"Dw1");
  cmd.mapvar(h,"h");

  de.yiGet(0,0)=w0;
  de.yiGet(0,1)=Dw0;
  de.yiGet(1,0)=w1;
  de.yiGet(1,1)=Dw1;
  de.hstep=h;

  double k0=0.1;
  double k1=0.03;
  double k2=0.2;

  cmd.mapvar(k0,"k0");
  cmd.mapvar(k1,"k1");
  cmd.mapvar(k2,"k2");
  de.integrator.yDorder.ki[0] = k0;
  de.integrator.yDorder.ki[1] = k1;
  de.integrator.yDorder.ki[2] = k2;

  string filename("myfile.txt");
  cmd.mapvar(filename,"filename");

  ofstream fil(filename.c_str(),ios::trunc);

  de.xval = 0.0;

  cout << de << endl;

  double xmax(0.3);
  cmd.mapvar(xmax,"xmax");

  uint width(14);
/*

  fil << setw(width) << de.xval;
  fil << setw(width) << de.yiGet(0,0);
  fil << setw(width) << de.yiGet(1,0);
  fil << endl;
*/

  long unsigned int sample=1;
  long unsigned int counter=0;

  cmd.mapvar(sample,"sample");

  for (; de.xval < xmax; de.eval(), ++counter )
  {
    if ((counter % sample)!=0)
      continue;
    

    fil << setw(width) << de.xval;
    fil << setw(width) << de.yiGet(0,0);
    fil << setw(width) << de.yiGet(1,0);
    fil << endl;
  }
//  cout << de << endl;


}