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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
|
#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 */
/*cout<<"enter the initial value for variable as (ti)"<<endl;
cin>>ti;
cout<<"initial value for function x(t) as xi"<<endl;
cin>>xi;
cout<<"enter the initial (vi)"<<endl;
cin>>vi;
cout<<"enter the step size for integration as dt"<<endl;
cin>>dt;
cout<<"enter the tmax to integrate from ti till tmax"<<endl;
cin>>tmax; */
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;
}
|
