1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
|
//this code works well//
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;
/* function prototypes */
double f1(double, double, double);
double f2(double, double, double);
double euler2m(double(*)(double, double, double),
double(*)(double, double, double),
double, double, double, double,
double&, double&);
int main()
{
double ti, xi, vi, tf, xf, vf, dt, tmax;
/* initial information */
ifstream inputFile1("afile-ti.txt");
ifstream inputFile2("afile-xi.txt");
ifstream inputFile3("afile-vi.txt");
ifstream inputFile4("afile-dt.txt");
ifstream inputFile5("afile-tmax.txt");
ofstream outputFile("output3.txt");
/*
ti = 0.0; // initial value for variable
xi = 0.0; // initial value for function x(t)
vi = 1.0; // initial
dt = 0.1; // step size for integration
tmax = 12.0; // integrate from ti till tmax */
inputFile1>>tf;
inputFile2>>xi;
inputFile3>>vi;
inputFile4>>dt;
inputFile5>>tmax;
/* integration of ODE */
while (ti <= tmax)
{
tf = ti + dt;
euler2m(f1,f2,ti,xi,vi,tf,xf,vf);
ti = tf;
xi = xf;
vi = vf;
}
cout<<"t="<<tf<<"\t"<<"x="<<xf<<"\t"<<"x'="<<vf<<"\n";
outputFile<<"t="<<tf<<"\t"<<"x="<<xf<<"\t"<<"x'="<<vf<<"\n";
return 0;
}
/*
Definition of the x'(t) = f1(t,x,x') = x' by the definition
*/
double f1(double t, double x, double v)
{
double d1x;
d1x = v;
return d1x;
}
/*
* Definition of the x"(t) = f2(t,x,x')
*/
double f2(double t, double x, double v)
{
double d2x;
d2x = (t*t);
return d2x;
}
//-----------------------------------------------------------------------
double euler2m(double(*d1x)(double, double, double),
double(*d2x)(double, double, double),
double ti, double xi, double vi, double tf,
double& xf, double& vf)
{
xf = xi + d1x(ti,xi,vi)*(tf-ti);
vf = vi + d2x(ti,xi,vi)*(tf-ti);
/* correction */
xf = xi + (d1x(ti,xi,vi)+d1x(ti,xf,vf))*0.5*(tf-ti);
vf = vi + (d2x(ti,xi,vi)+d2x(ti,xf,vf))*0.5*(tf-ti);
return 0.0;
}
|