/*
  In this example we solve ordinary differential equation by 
  program dumka3.cpp (see http://www.math.tulane.edu/~amedovik
  y'=-500*(y-cos(t)), y(0)=0  
  This is stiff ODE and expected smooth solution should be approximately
  y-cos(t)=0
  Indeed, exact solution is
  y=-l^2/(l^2+1.)*exp(-l*t)+(l^2*cos(t)-l*sin(t))/(l^2+1), where l=500
  which is approximately cos(t)

  class VectorFunction - should have method 
                         Vector compute(Vector y, double t, Vector& result)
 
  isStepAccepted - is not necesseraly in general, but I've used it to
  pass information into the class that step has been accepted, because 
  dumka3 rejects step automatically if error is larger than tolerance

  cour - template-function, which returns lower evaluation of maximal step size 
         of explicit Euler method
         for example in equation y'=-500(y-cos(t))
         maximal step size is 2/500, so I can take 2./500 or
         something little bit smaller like 2./500.1
*/
//Command for SGI:
//CC oneStiffODE.cpp -lm -LANG:std -64 -apo

#include <stdio.h>
#include <math.h>
#include <iostream>
#include <fstream>
//#include "MatrixVectorVectorFunction.h" // you can use my Vector and Matrix classes
#include <vector>
#include "dumka3.cpp"
//#include "euler.cpp"


using namespace std;
ofstream  solution( "solution" ); // open output file for solution
ofstream  error( "error" ); // open output file for error
ofstream  exact( "exact" ); // open output file for exact solutions

template<class Vector, class Element>
class VectorFunction: public VectorFunctionBase<Vector,Element>{
  public:  
    bool isStepAccepted;
    VectorFunction() {
      isStepAccepted=false; 
    };
    
    // compute - calculates right hand side of the equations dy/dt=f(y,t),and 
    //           returns "f(y,t)" in "result" argument
    void compute (Vector& y, double t, Vector& result) {
      result[0]= -500.*(y[0]-cos(t));
      if(isStepAccepted) {
        isStepAccepted=false;
        solution << t << " " << y[0] << endl;
        double _exact = (-(500.*500.)/(500.*500.+1.))*exp(-500.*t)+
              ((500.)*(500.)*cos(t)-500.*sin(t))/(500.*(500.)+1.);
        error << t  << " " << _exact - y[0] << endl;
        exact << t << " " << _exact << endl; 
      };
    };

  double cour(Vector& y, double t) {
    double cou=2./500.1;
    return cou;
  };
  Element GetNonZeroElement() {
    return rand()%100/100.;
  }
};



int main(int argc  , char** argv ){
  vector<double>  y(1);
  VectorFunction<vector<double>, double>f;
  y[0]=0.;            // y - solution of ODE
  double tinitial=0.;
  double t=tinitial;    // t - initial time
  double tend=1.;       // tend - interval of integration
  double rtol=0.1; // rtol - relative error 
  double atol=0.1;      // tolerance of computation
  double h0=0.1*f.cour(y,t);     // h0 - initial step size
                      // f - class VectorFunction which has 
                      //     method Vector compute(Vector , double )
  dumka3(t, tend, h0, atol, rtol, f, y, true );
  //euler( t, tend, f, y );
  solution.close();   // close output file
  error.close();
};