ode/ode1.cpp

   1 #include <cmath>
   2 #include <cstdlib>
   3 #include <iostream>
   4 #include <vector>
   5
   6 using namespace std;
   7
   8 class odeSolver {
   9 public:
  10     odeSolver() {}
  11     void addY(vector<float> *y) { this->y = y; }
  12     void addDYDX(vector<float> *dydx) { this->dydx = dydx; }
  13     void setDerivs(void (*pderivs)(float, vector<float> *, vector<float> *));
  14     void printY(void);
  15     virtual void getNextCalc(void) {};
  16     void setN(unsigned int n) { this->n = n; }
  17     void setStep(float step) { this->step = step; }
  18     void setX(float x) { this->x = x; }
  19 protected:
  20     float step;
  21     float x;
  22     void (*pderivs)(float, vector<float> *, vector<float> *);
  23     vector<float> *y;
  24     vector<float> *dydx;
  25     unsigned int n;
  26 };
  27
  28 void odeSolver::printY(void)
  29 {
  30     vector<float>::iterator it;
  31
  32     for (it = y->begin(); it != y->end(); it++)
  33         cout << *it << endl;
  34 }
  35 void odeSolver::setDerivs(void (*pderivs)(float,
  36                                vector<float> *,
  37                                vector<float> *))
  38 {
  39     this->pderivs = pderivs;
  40 }
  41
  42 // Based on code from "Numerical Recipes in C, 2nd Ed."
  43 class odeSolverRK4 : public odeSolver
  44 {
  45 public:
  46     void getNextCalc(void)
  47     {
  48         unsigned int i;
  49         float xh, hh, h6;
  50         vector<float> dym(n), dyt(n), yt(n);
  51
  52         hh = step * 0.5;
  53         h6 = step / 6.0;
  54         xh = x + hh;
  55
  56         //
  57         for (i = 0; i < n; i++)
  58             yt[i] = (*y)[i] + hh * (*dydx)[i];
  59         pderivs(xh, &yt, &dyt);
  60
  61         //
  62         for (i = 0; i < n; i++)
  63             yt[i] = (*y)[i] + hh * dyt[i];
  64         pderivs(xh, &yt, &dym);
  65
  66         //
  67         for (i = 0; i < n; i++) {
  68             yt[i] = (*y)[i] + step * dym[i];
  69             dym[i] += dyt[i];
  70         }
  71         pderivs(x + step, &yt, &dyt);
  72
  73         //
  74         for (i = 0; i < n; i++)
  75             (*y)[i] = (*y)[i] + h6 * ((*dydx)[i] + dyt[i] + 2.0 * dym[i]);
  76
  77         //
  78         x += step;
  79     }
  80 };
  81
  82 void myderivs(float x, vector<float> *y, vector<float> *dydx)
  83 {
  84
  85     (*dydx)[0] = x * cos(x * (*y)[0]);
  86 }
  87
  88 int main(void)
  89 {
  90     odeSolverRK4 myos;
  91     vector<float> y;
  92     vector<float> dydx;
  93
  94     // y vector
  95     y.push_back(1.0);
  96     myos.addY(&y);
  97
  98     // dydx vector
  99     dydx.push_back(1.0);
 100     myos.addDYDX(&dydx);
 101
 102     myos.setDerivs(myderivs);
 103     myos.setN(1);
 104     myos.setX(0);
 105     myos.setStep(0.07);
 106
 107     for (int i = 0; i < 1000; i++) {
 108         myos.getNextCalc();
 109         myos.printY();
 110     }
 111
 112     return EXIT_SUCCESS;
 113 }