数值代数–理查森外推法
实验四
一、实验名称
理查森外推算法
二、实验目的与要求:
实验目的:掌握理查森外推算法。
实验要求:1. 给出理查森外推算法思路,
2. 用C语言实现算法,运行环境为Microsoft Visual C++。
三、算法思路:
1. 假设函数泰勒展开式可表示为 和,将两式相减,消去偶数项,则,整理得到下式,记L表示,表示微分形式,则有 (1)用h/2代替h,有 (2),由(1)(2)两式子有推广这种方法,就是理查森外推法了。
2. 理查森外推法公式 , , 用下列公式计算,k=1,2,…,M,n=k,k+1,…,M。 则有,当n和k足够大时D(n,k)可充分接近。
3. 上机算法 input h , M
for n=0 to M do
D(n , 0)
end do
for k=1 to M do
for n=k to M do end do
end do
output D(n , k)
四、实验题目:
五、问题的解:
编写程序(程序见后面附录),输出结果如下:
分析得到的结果,发现在对角线附近D(n , k)的值越来越稳定,通过上面算法阐述,我们知道D(n , k)应该是越来越接近我们想求到的导数的,与实验结果一致。
六、附录:
实验编程,运行环境为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");
}
}