nintegration.h

Go to the documentation of this file.

#ifndef NINTEGRATION_H
#define NINTEGRATION_H

#include <cassert>
#include <iostream>
using namespace std;

#include <print.h>

template< typename D, typename T >
class nintegrationEuler 
{
public:

  D yDorder;

  nintegrationEuler(T * xi_, T * yi_)
    : yDorder(xi_,yi_) {}

//  void operator ++ ()
//    { ++yDorder; }

  void eval
  ( 
    T & dz,
    T const h,
    T const z, 
    T const x 
  )
    { yDorder(dz,z,x); dz*=h; }
};


template< typename D, typename T >
class nintegrationTrapezoidEuler
{
public:

  D yDorder;

  nintegrationTrapezoidEuler(T * xi_, T * yi_)
    : yDorder(xi_,yi_) {}

  void operator ++ ()
    { ++yDorder; }

  void eval
  ( 
    T & dy,
    T const h,
    T const y, 
    T const x 
  )
  {
    // Estimate the implicit variable y[n+1] by Eulers method.
    // w1 is y[n+1] interpolated: y[n+1] = y[n]+h*y'[n]

    T yD;
    yDorder(yD,y,x);
    T w1 = y + h*yD;

    T x1 = x+h;
    T y1D;
    yDorder(y1D,w1,x1);

    dy = (y1D+yD)*h/2.0;
  }

};


template< typename D, typename T >
class nintegrationRK
{
public:

  D yDorder;

  nintegrationRK(T * xi_, T * yi_)
    : yDorder(xi_,yi_) {}

  void operator ++ ()
    { ++yDorder; }

  void eval
  ( 
    T & dy,
    T const h,
    T const y, 
    T const x 
  )
  {
    T c1;
    yDorder(c1,y,x);
    c1 *= h;

    T x1 = x + h*0.5;
    T y1 = y + c1*0.5;

    T c2;
    yDorder(c2,y1,x1);
    c2 *= h;

    T y2 = y+c2*0.5;
    T c3;
    yDorder(c3,y2,x1);
    c3 *= h;

    T x2 = x + h;
    T y3 = y + c3;
    T c4;
    yDorder(c4,y3,x2);
    c4 *= h;

    dy = (c1+c2*2.0+c3*2.0+c4)/6.0;
  }

};

#endif