实验四一、实验名称理查森外推算法二、实验目的与要求 ：实验目的：掌握理查森外推算法。实验要求： 1. 给出理查森外推算法思路，2. 用 C语言实现算法， 运行环境为 Microsoft Visual C++。三、算法思路：1. 假设函数)(xf泰勒展开式可表示为0)()(!1)(kkkxfhkhxf和0)()()1(!1)(kkkkxfhkhxf，将两式相减， 消去偶数项， 则0)12(12)()!12(2)()(kkkxfhkhxfhxf，整理得到下式1)12(2')()!12(1)]()([21)(kkkxfhkhxfhxfhxf，记 L表示)(' xf，)(h表示微分形式)]()([21hxfhxfh，则有4422)(Lhahah（1）用 h/2 代替 h，有16/4/)2/(L4422hahah（2），由（1）（2）两式子有16/4/)(31)2/(34L6644hahahh推广这种方法，就是理查森外推法了。2. 理查森外推法公式)2/()0,(nhnD, Mn0, 用下列公式计算)1,1(141)1,(144),(knDknDknDkkk，k=1,2,⋯,M，n=k,k+1,⋯，M 。 则有)(),(2khOLknD，当 n 和 k 足够大时 D(n,k)可充分接近)(' xf。3. 上机算法input h , M for n=0 to M do D(n , 0))2/(nhend do for k=1 to M do for n=k to M do )14/()]1,1()1,([)1,(),(kknDknDknDknDend do end do output D(n , k) )0,0(nkMn四、实验题目：五、问题的解：编写程序（程序见后面附录） ，输出结果如下：分析得到的结果，发现在对角线附近D(n , k)的值越来越稳定，通过上面算法阐述，我们知道D(n , k)应该是越来越接近我们想求到的导数)(' xf的，与实验结果一致。六、附录：实验编程，运行环境为Microsoft Visual C++ #include #include #include double f1(double x) //定义函数 f1(x)// { double y; y=(log(3.0+x)-log(3.0-x))/(2.0*x); return(y); } double f2(double x) //定义函数 f2(x)// { double y; y=(tan(asin(0.8)+x)-tan(asin(0.8)-x))/(2.0*x); return(y); } double f3(double x) //定义函数 f3(x)// { double y; y=(sin(x*x+x/3.0)-sin(x*x-x/3.0))/(2.0*x); return(y); } void main() { double D1[4][4],D2[5][5],D3[6][6]; int i,j; for(i=0;i<=3;i++) /* 第一个问题的理查森算法*/ D1[i][0]=f1(1.0/pow(2,i)); for(j=1;j<=3;j++) for(i=j;i<=3;i++) D1[i][j]=D1[i][j-1]+(D1[i][j-1]-D1[i-1][j-1])/(pow(4,j)-1); printf(" 第一道题结果： \n"); for(i=0;i<=3;i++) {for(j=0;j<=i;j++) printf("%0.12f ",D1[i][j]); printf("\n"); } for(i=0;i<=4;i++) /* 第二个问题的理查森算法*/ D2[i][0]=f2(1.0/pow(2,i)); for(j=1;j<=4;j++) for(i=j;i<=4;i++) D2[i][j]=D2[i][j-1]+(D2[i][j-1]-D2[i-1][j-1])/(pow(4,j)-1); printf(" 第二道题结果： \n"); for(i=0;i<=4;i++) {for(j=0;j<=i;j++) printf("%0.12f ",D2[i][j]); printf("\n"); } for(i=0;i<=5;i++) /* 第三个问题的理查森算法*/ D3[i][0]=f3(1.0/pow(2,i)); for(j=1;j<=5;j++) for(i=j;i<=5;i++) D3[i][j]=D3[i][j-1]+(D3[i][j-1]-D3[i-1][j-1])/(pow(4,j)-1); printf(" 第三道题结果： \n"); for(i=0;i<=5;i++) {for(j=0;j<=i;j++) printf("%0.12f ",D3[i][j]); printf("\n"); } }
