tmain.cpp - numeric - C++ library with numerical algorithms
 (HTM) git clone git://src.adamsgaard.dk/numeric
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
       tmain.cpp (1370B)
       ---
            1 #include <iomanip>
            2 #include <cstdlib> // Exit function
            3 #include <string.h>
            4 #include <cmath>
            5 #include <vector>
            6 
            7 // File I/O
            8 #include <iostream> 
            9 #include <fstream>
           10 
           11 using namespace std;
           12 
           13 // this function gives the sign
           14 bool sign(double n)
           15 { return n >= 0; };
           16 
           17 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 );
           18 
           19 // this function simulates a cannonball
           20 void func1(int n, double x, vector<double>* y, vector<double>* dydx){
           21   int A = sign((*y)[2]);
           22   int B = sign((*y)[3]);
           23   double k = 0.02;
           24   double g = 9.82;
           25   (*dydx)[0]=(*y)[2];
           26   (*dydx)[1]=(*y)[3];
           27   (*dydx)[2]=sqrt((*y)[2]*(*y)[2]+(*y)[3]*(*y)[3])*-A*k;
           28   (*dydx)[3]=-g-sqrt((*y)[2]*(*y)[2]+(*y)[3]*(*y)[3])*B*k;
           29 };
           30 
           31 int main () {
           32 
           33   int n=4;  // Number of equations
           34   int max=10000;  // Maximum "iterations"
           35   double a=0, b=15, h=0.01, acc=0.001, eps=0.001;
           36   vector<double> tlist(max);  // 
           37   vector<vector<double> > ylist(max,vector<double>(n)); 
           38   tlist[0]=a; ylist[0][0]=0; ylist[0][1]=0; ylist[0][2]=100;
           39   ylist[0][3]=100;
           40   
           41   rkdriver(&func1,n,&tlist,&ylist,b,h,acc,eps,max); 
           42   
           43   // Printing the output
           44   int it = 0;
           45   while(tlist[it]<b) {
           46     cout << tlist[it] << " ";
           47     for (int j = 0; j<n;j++) cout << ylist[it][j] << " ";
           48     cout << endl;
           49     it++;
           50   };        
           51   
           52 };
           53