数值代数-理查森外推法_百度文库.DOCVIP

实验四 一、实验名称 理查森外推算法 二、实验目的与要求： 实验目的：掌握理查森外推算法。 实验要求：1. 给出理查森外推算法思路， 2. 用C语言实现算法，运行环境为Microsoft Visual C++。 三、算法思路： 1 1. 假设函数f(x)泰勒展开式可表示为f(x?h)??hkf(k)(x) k?0?k! 1kk(k)f(x?h)?(?1)hf(x)，将两式相减，消去偶数项，则和? k?0k!? f(x?h)?f(x?h)?? 2h2k?1f(2k?1)(x)，整理得到下式k?0(2k?1)!??11f(x)?[f(x?h)?f(x?h)]??h2kf(2k?1)(x)，记L2hk?1(2k?1)! 1表示f(x)，?(h)表示微分形式[f(x?h)?f(x?h)]，则2h 24L??(h)?ah?ah?? （1）用h/2代替h，有有24 L??(h/2)?a2h2/4?a4h4/16?? （2），由（1）（2）两式4146L??(h/2)??(h)?ah/4?ah/16??推广这种46子有33 方法，就是理查森外推法了。 2. 理查森外推法公式D(n,0)??(h/2n) , 0?n?M, 用下4k1D(n,k?1)?kD(n?1,k?1)，列公式计算D(n,k)?k4?14?1 2kD(n,k)?L?O(h)，k=1,2,…,M，n=k,k+1,…，M。 则有 f当n和k足够大时D(n,k)可充分接近(x)。 3. 上机算法 input h , M for n=0 to M do n?(h/2) D(n , 0) ? end do for k=1 to M do for n=k to M do D(n,k)?D(n,k?1)?[D(n,k?1)?D(n?1,k?1)]/(4k?1) end do end do output D(n , k) (0?n?M,0?k?n) 四、实验题目： 五、问题的解： 编写程序（程序见后面附录），输出结果如下： 分析得到的结果，发现在对角线附近D(n , k)的值越来越稳定，通过上面算法阐述，我们知道D(n , k)应该是越来越接近我们想求到的导数 六、附录： 实验编程，运行环境为Microsoft Visual C++ #include math.h #include stdio.h #include stdlib.h double f1(double x) //定义函数f1(x)// { double y; y=(log(3.0+x)-log(3.0-x))/(2.0*x); return(y); } f(x)的，与实验结果一致。 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[