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
|
/
#define N 2
//#define dist 0.1
#define MAX 30.0
FILE *output;
void runge4(double x, double y[], double step);
double f(double x, double y[], int i);
main()
{
double t, y[N];
double dist;
output=fopen("osc.dat", "w");
y[0]=1.0;
y[1]=0.0;
fprintf(output, "0\t%f\n", y[0]);
int j=1;
double s=0;
while (dist*j<=MAX)
{
t= j*dist;
runge4(t,y,dist);
j++;
if (y[0]<0)
break;
fprintf(output, "%f\t%f\n", t, y[0]);
}
fclose(output);
}
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)
{
............
}
|