trkdriver.cpp - numeric - C++ library with numerical algorithms
 (HTM) git clone git://src.adamsgaard.dk/numeric
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
       trkdriver.cpp (1286B)
       ---
            1 #include <math.h>
            2 #include <stdlib.h>
            3 #include <stdio.h>
            4 #include <iostream> 
            5 #include <fstream>
            6 #include <vector>
            7 using namespace std;
            8 
            9 void rkstep(void f(int,double,vector<double>*,vector<double>*),int n, double t,vector<double>* y,double h,vector<double>* y1,vector<double>* dy);
           10 
           11 // Main driver function
           12 void rkdriver(void f(int, double, vector<double>*, vector<double>*),int n, vector<double>* tlist, vector<vector<double> >* ylist, double b, double h, double acc, double eps, int max ) {
           13 
           14   int i=0;  //iterator
           15   double t = (*tlist)[0];
           16   double a = t; 
           17   vector<double> dy(n);
           18   vector<double> y1(n);
           19   
           20   while((*tlist)[i]<b) {
           21     double t=(*tlist)[i];
           22     vector<double> y=(*ylist)[i];
           23     if(t+h>b) h=b-t;
           24     rkstep(f,n,t,&y,h,&y1,&dy);
           25     double err=0; for(int j=0;j<n;j++)err+=dy[j]*dy[j];
           26     err=sqrt(err);
           27     double normy=0; for(int j=0;j<n;j++)normy+=y1[j]*y1[j];
           28     normy=sqrt(normy);
           29     double tol=(normy*eps+acc)*sqrt(h/(b-a));
           30 
           31     if(tol>err){ // accept step and go on
           32       i++;
           33       if(i>max-1) { 
           34         cout << "Reached max step \nIncrease max step to go further \n";
           35         exit(1);
           36       };
           37 
           38       (*tlist)[i]=t+h;
           39       
           40       for(int j=0;j<n;j++) (*ylist)[i][j]=y1[j];
           41     };
           42 
           43     if(err>0) h = h*pow(tol/err,0.25)*0.95;
           44     else h = 2*h;
           45   }//end while                                
           46 };