GopherProxy

	trkdriver.cpp - numeric - C++ library with numerical algorithms
	git clone git://src.adamsgaard.dk/numeric
	Log
	Files
	Refs
	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, do…
	10
	11 // Main driver function
	12 void rkdriver(void f(int, double, vector<double>, vector<double>),int …
	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=(normyeps+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 = hpow(tol/err,0.25)0.95;
	44 else h = 2*h;
	45 }//end while
	46 };