#include <stdio.h>
#include <math.h>
#define N 2 /* number of first order equations */
#define dist 0.1 /* stepsize in t*/
#define MAX 30.0 /* max for t */
FILE
*output;
FILE
*output1;
FILE
*pipe;
void
runge4(
double
x,
double
y[],
double
step);
double
f(
double
x,
double
y[],
int
i);
void
main(){
double
t, y[N];
int
j;
output=
fopen
(
"osc.dat"
,
"w"
);
output1=
fopen
(
"osc1.dat"
,
"w"
);
y[0]=1.0;
y[1]=0.0;
for
(j=1; j*dist<=MAX ;j++)
{
t=j*dist;
runge4(t, y, dist);
fprintf
(output,
"%f\t%f\n"
, t, y[0]);
fprintf
(output1,
"%f\t%f\n"
, t, y[1]);
}
fclose
(output);
fclose
(output1);
pipe = popen(
"gnuplot -persist"
,
"w"
);
fprintf
(pipe,
"set term png enhanced font \"y:/wqy-microhei.ttc\" 18 \n"
);
fprintf
(pipe,
"set output \"test.png\"\n"
);
fprintf
(pipe,
"plot \"osc.dat\" title \"位移\" with lines, \"osc1.dat\" title \"速度\" with lines\n"
);
fprintf
(pipe,
"quit\n"
);
fprintf
(pipe,
"quit\n"
);
pclose(pipe);
}
void
runge4(
double
x,
double
y[],
double
step){
double
h=step/2.0,
t1[N], t2[N], t3[N],
k1[N], k2[N], k3[N],k4[N];
int
i;
for
(i=0;i<N;i++){
t1[i]=y[i]+0.5*(k1[i]=step*f(x,y,i));
}
for
(i=0;i<N;i++){
t2[i]=y[i]+0.5*(k2[i]=step*f(x+h, t1, i));
}
for
(i=0;i<N;i++){
t3[i]=y[i]+ (k3[i]=step*f(x+h, t2, i));
}
for
(i=0;i<N;i++){
k4[i]= step*f(x+step, t3, i);
}
for
(i=0;i<N;i++){
y[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
}
}
double
f(
double
x,
double
y[],
int
i){
if
(i==0)
x=y[1];
if
(i==1)
x=-y[0]-0.5*y[1];
return
x;
}