数值代数-理查森外推法

实验四一、实验名称理查森外推算法二、实验目的与要求：实验目的：掌握理查森外推算法。实验要求：1.给出理查森外推算法思路，2.用C语言实现算法，运行环境为MicrosoftVisualC++。三、算法思路：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,Mn0,用下列公式计算)1,1(141)1,(144),(knDknDknDkkk，k=1,2,…,M，n=k,k+1,…，M。则有)(),(2khOLknD，当n和k足够大时D(n,k)可充分接近)('xf。3.上机算法inputh,Mforn=0toMdoD(n,0))2/(nhenddofork=1toMdoforn=ktoMdo)14/()]1,1()1,([)1,(),(kknDknDknDknDenddoenddooutputD(n,k))0,0(nkMn四、实验题目：五、问题的解：编写程序（程序见后面附录），输出结果如下：分析得到的结果，发现在对角线附近D(n,k)的值越来越稳定，通过上面算法阐述，我们知道D(n,k)应该是越来越接近我们想求到的导数)('xf的，与实验结果一致。六、附录：实验编程，运行环境为MicrosoftVisualC++#includemath.h#includestdio.h#includestdlib.hdoublef1(doublex)//定义函数f1(x)//{doubley;y=(log(3.0+x)-log(3.0-x))/(2.0*x);return(y);}doublef2(doublex)//定义函数f2(x)//{doubley;y=(tan(asin(0.8)+x)-tan(asin(0.8)-x))/(2.0*x);return(y);}doublef3(doublex)//定义函数f3(x)//{doubley;y=(sin(x*x+x/3.0)-sin(x*x-x/3.0))/(2.0*x);return(y);}voidmain(){doubleD1[4][4],D2[5][5],D3[6][6];inti,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);}}