c语言差分求导,急求:如何解决高阶微分求导的计算误差问题！！！

该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

需寻求解高阶微分(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;

}