double f1 double double double double step double gamma double omega r

 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
double f1 (double t, double x, double y,double step,double gamma,double omega )
{
return step * ( -2 * gamma * y - pow(omega,2) * x) ;
}
double f2 (double t, double x, double y,double step,double gamma,double omega )
{
return step * y ;
}
int main() {
const double time =100;
const double step =0.01;
const unsigned int count = int((time)/step) + 1;
vector<double> time_vec(count);
//vector<double> new_x(count);
//vector<double> new_y(count);
double X[(int)count];
double X1[(int)count];
double X2[(int)count];
double X3[(int)count];
double X4[(int)count];
double Y1[(int)count];
double Y2[(int)count];
double Y3[(int)count];
double Y4[(int)count];
double Y[(int)count];
double t[(int)count];
X[0]=5;Y[0] =5;t[0]=0;
for(int i= 1 ; i < count ; i++)
{
double ntime = i * step ;
X1[i]=f2(t[i-1],X[i-1],Y[i-1], step, gamma, omega);
Y1[i]=f1(t[i-1],X[i-1],Y[i-1], step, gamma, omega);
X2[i]=f2(t[i-1] + step/2, X[i-1]+ X1[i]/2.0,Y[i-1]+Y1[i]/2.0 ,step, gamma, omega);
Y2[i]=f1(t[i-1] + step/2, X[i-1]+ X1[i]/2.0,Y[i-1]+Y1[i]/2.0 ,step, gamma, omega);
X3[i]=f2(t[i-1] + step/2, X[i-1]+ X2[i]/2,Y[i-1]+Y2[i]/2,step, gamma, omega);
Y3[i]=f1(t[i-1] + step/2, X[i-1]+X2[i]/2,Y[i-1]+Y2[i]/2,step, gamma, omega);
X4[i]=f2(t[i-1] + step/2, X[i-1]+X3[i],Y[i-1]+Y3[i],step, gamma, omega);
Y4[i]=f1(t[i-1] + step/2, X[i-1]+X3[i],Y[i-1]+Y3[i],step, gamma, omega);
X[i]=X[i-1]+(X1[i]+2*X2[i]+2*X3[i]+X4[i])/6;
Y[i]=Y[i-1]+(Y1[i]+2*Y2[i]+2*Y3[i]+Y4[i])/6;
time_vec[i] = ntime;
}