该楼层疑似违规已被系统折叠 隐藏此楼查看此楼
需寻求解高阶微分(5阶)的程序,我试着用中心差分法做了个简单的,
但是如果直接求已知函数(如sin(x)),除了边界有些误差外,还马马虎虎,
但是如果用读取的函数值文件(如“x_y(N=4096)-1.data”)的话,从
2阶导数开始就出现严重误差了,请教:
1,算法是否有误?
2,N的取值多少合适?
拜求各位高手了!!!!!!
-------------------------------------------------
#include
#include
/*-------------------------------------------------------*/
/*Differential (Central Difference Method) */
/*-------------------------------------------------------*/
void main()
{
#define N(4096)//2^12
longi;
double x[N],y[N],y1[N],y2[N],y3[N],y4[N],y5[N];
double delta;
double PI=acos(-1.0);
FILE *fh,*f0,*fh0;
if(N!=0)
{
delta=2.0*PI/(double)N;
#if 0
for (i=0; i
x[i]=i*delta;
y[i]=sin(x[i]);
}
#endif
#if 1
f0=fopen("x_y(N=4096)-1.data","r");
for (i= 0; i
fscanf(f0,"%lf %lf\n", &x[i], &y[i]);
}
fclose(f0);
#endif
for (i=1; i
{
y1[i]=(y[i+1]-y[i-1])/(2.0*delta);
}
y1[0]=(y[1]-y[0])/delta;
y1[N-1]=(y[N-1]-y[N-2])/delta;
for (i=1; i
{
y2[i]=(y1[i+1]-y1[i-1])/(2.0*delta);
}
y2[0]=(y1[1]-y1[0])/delta;
y2[N-1]=(y1[N-1]-y1[N-2])/delta;
for (i=1; i
{
y3[i]=(y2[i+1]-y2[i-1])/(2.0*delta);
}
y3[0]=(y2[1]-y2[0])/delta;
y3[N-1]=(y2[N-1]-y2[N-2])/delta;
for (i=1; i
{
y4[i]=(y3[i+1]-y3[i-1])/(2.0*delta);
}
y4[0]=(y3[1]-y3[0])/delta;
y4[N-1]=(y3[N-1]-y3[N-2])/delta;
for (i=1; i
{
y5[i]=(y4[i+1]-y4[i-1])/(2.0*delta);
}
y5[0]=(y4[1]-y4[0])/delta;
y5[N-1]=(y4[N-1]-y4[N-2])/delta;
}
//fh0=fopen("x_y(N=4096)-1.data","a");
fh=fopen("x_y_y1-5.data","a");
for (i=0; i
//fprintf(fh0,"%g %g\n",x[i],y[i]);
fprintf(fh,"%g %g %g %g %g %g %g\n",x[i],y[i],y1[i],y2[i],y3[i],y4[i],y5[i]);
}
//fclose(fh0);
fclose(fh);
return;
}