/* In this example we solve ordinary differential equation by program dumka3.cpp (see http://www.math.tulane.edu/~amedovik y'=-500*(y-cos(t)), y(0)=0 This is stiff ODE and expected smooth solution should be approximately y-cos(t)=0 Indeed, exact solution is y=-l^2/(l^2+1.)*exp(-l*t)+(l^2*cos(t)-l*sin(t))/(l^2+1), where l=500 which is approximately cos(t) class VectorFunction - should have method Vector compute(Vector y, double t, Vector& result) isStepAccepted - is not necesseraly in general, but I've used it to pass information into the class that step has been accepted, because dumka3 rejects step automatically if error is larger than tolerance cour - template-function, which returns lower evaluation of maximal step size of explicit Euler method for example in equation y'=-500(y-cos(t)) maximal step size is 2/500, so I can take 2./500 or something little bit smaller like 2./500.1 */ //Command for SGI: //CC oneStiffODE.cpp -lm -LANG:std -64 -apo #include #include #include #include //#include "MatrixVectorVectorFunction.h" // you can use my Vector and Matrix classes #include #include "dumka3.cpp" //#include "euler.cpp" using namespace std; ofstream solution( "solution" ); // open output file for solution ofstream error( "error" ); // open output file for error ofstream exact( "exact" ); // open output file for exact solutions template class VectorFunction: public VectorFunctionBase{ public: bool isStepAccepted; VectorFunction() { isStepAccepted=false; }; // compute - calculates right hand side of the equations dy/dt=f(y,t),and // returns "f(y,t)" in "result" argument void compute (Vector& y, double t, Vector& result) { result[0]= -500.*(y[0]-cos(t)); if(isStepAccepted) { isStepAccepted=false; solution << t << " " << y[0] << endl; double _exact = (-(500.*500.)/(500.*500.+1.))*exp(-500.*t)+ ((500.)*(500.)*cos(t)-500.*sin(t))/(500.*(500.)+1.); error << t << " " << _exact - y[0] << endl; exact << t << " " << _exact << endl; }; }; double cour(Vector& y, double t) { double cou=2./500.1; return cou; }; Element GetNonZeroElement() { return rand()%100/100.; } }; int main(int argc , char** argv ){ vector y(1); VectorFunction, double>f; y[0]=0.; // y - solution of ODE double tinitial=0.; double t=tinitial; // t - initial time double tend=1.; // tend - interval of integration double rtol=0.1; // rtol - relative error double atol=0.1; // tolerance of computation double h0=0.1*f.cour(y,t); // h0 - initial step size // f - class VectorFunction which has // method Vector compute(Vector , double ) dumka3(t, tend, h0, atol, rtol, f, y, true ); //euler( t, tend, f, y ); solution.close(); // close output file error.close(); };