desys.h

Go to the documentation of this file.

#ifndef DESYS_H
#define DESYS_H

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

#include <print.h>

template< typename NINT, typename T >
class desys
{
  void cleanup();
  desys()
    { assert(false); }
public:

  uint dim;
  uint order;

  T * dy;

  T xval;
  T * yi;

  T & yiGet(uintc i, uintc k)
    { assert(i<dim); assert(k<order); return yi[k+i*order]; }

  T const yiGet(uintc i, uintc k) const
    { assert(i<dim); assert(k<order); return yi[k+i*order]; }

  NINT integrator;

  T hstep;

  desys(uintc dim_, uintc order_);

  ~desys();

  void eval();

  ostream & print(ostream & os) const;
};


template< typename NINT, typename T >
ostream & operator << (ostream & os, desys<NINT,T> const & x)
  { return x.print(os); }

//---------------------------------------------------
//  Implementation


template< typename NINT, typename T >
ostream & desys<NINT,T>::print(ostream & os) const
{
  os << "dim=" << dim << endl;
  os << "order=" << order << endl;
  os << "hstep=" << hstep << endl;
  os << "xval=" << xval << endl;
  os << "(y0, y0', ...)" << endl;
  os << "(y1, y1', ...)" << endl;
  os << "=" << endl;

  os << "yi" << endl;
  for (uint i=0; i<dim; ++i)
  {
    os << "[" << i << "]: ";
    for (uint k=0; k<order; ++k)
      os << yi[i*order+k] << " ";
    os << endl;
  }
 
  return os;
}


template< typename NINT, typename T >
desys<NINT,T>::~desys()
{
  delete[] dy;
  delete[] yi;
}

template< typename NINT, typename T >
desys<NINT,T>::desys(uintc dim_, uintc order_)
  : dim(dim_), order(order_), dy(new T[dim]),
    yi(new T[order*dim]), integrator(&xval,yi) 
{
  assert(dim>0);
  assert(order>0);
  xval=0.0;
}


template< typename NINT, typename T >
void desys<NINT,T>::eval()
{
  uint c;
  assert(hstep!=0);
  for (uint i=0; i<dim; ++i )
  {
    c = i*order+order-1;

//cout << SHOW(xval) << endl;

//cout << SHOW(order) << endl;
//cout << SHOW(c) << endl;
//cout << SHOW(order+c) << endl;
//cout << SHOW(i*(order+1)) << endl;
//cout << SHOW(yi[c]) << endl;
    integrator.eval(dy[i],hstep,yi[c],xval); 
//print(cout); cout << endl;
    yi[c] += dy[i];
//cout << SHOW(yi[c]) << endl;
//cout << SHOW(yi[order-1+c]) << endl;
//cout << "*";
//print(cout); cout << endl;

    ++integrator.yDorder;
  }
  xval += hstep;

  if (order<=1)
    return;
    
  for (uint i=0; i<dim; ++i )
  {
    c = i*order + order;

    for (uint k=2; k<=order; ++k ) 
      yi[c-k] += yi[c-k+1]*hstep;
  }
}








#endif