lichaxun.c

理查逊外推法求数值微分

#include<stdio.h>
#include<stdio.h>
#include<math.h>
#define MN 32
int RN=0;
double EPS=0.4E-12,ERR,ANS;
double H[MN],T0[MN],T1[MN],T2[MN],T3[MN];
double x0=4.0;
double h0=0.32;
double f(double x){return (x)*(x)*(x);}
double FT0(double x,double h)
{return (f(x+h)-f(x-h))/h/2;}
int Operation()
{
	int k;
	H[0]=h0;
	T0[0]=FT0(x0,h0);
	T1[0]=T0[0];
	T2[0]=T0[0];
	T3[0]=T0[0];
	for(k=1;k<MN;k++)
	{H[k]=H[k-1]/2.0;
	T0[k]=FT0(x0,H[k]);
	T1[k]=(4.0*T0[k]-T0[k-1])/3.0;
	T2[k]=(16.0*T1[k]-T1[k-1])/15.0;
	T3[k]=(64.0*T2[k]-T2[k-1])/63.0;
	if(k<4)continue;
	ERR=fabs(T3[k]-T3[k-1]);
	if(ERR<EPS)break;
	}
	if(k<MN)k++;
	RN=k;
	ANS=T3[RN-1];
	return RN;
}
   
   main()
   {Operation();
   system("cls");
   printf("\n   ANS=%12.8lf\n   ERR=%12.8lf\n",ANS,ERR);
   }