数值666 (2)
说明:文章内容仅供预览,部分内容可能不全。下载后的文档,内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的,是否完整无缺。
龙贝格求积分算法 4.龙贝格求积分算法 4.1算法说明 生成JK的逼近表R(J,K),并以R(J1,J1)为最终解来逼近积分 baf(x)dxR(J,J) 逼近R(J,K)存在于一个特别的下三角矩阵中,第0列元素R(J,0)用基于2J个[a,b]子区间的连续梯形方法计算,然后利用龙贝格公式计算R(J,K)。当1KJ时,第J行的元素为 R(J,K)R(J,K1)R(J,K1)R(J1,K1) 4K1当|R(J,J)R(J1,J1)|tol时,程序在第(J1)行结束。 4.2龙贝格求积分算法流程图 图4-1 流程图 第 14 页 共32页 数值计算课程设计 4.3龙贝格求积分算法程序调试 求解积分方程果如图4-2: 91xdx为例,对龙贝格求积分算法程序进行编译和链接,执行后得到结 4.4龙贝格求积分算法程序代码 图4-2 执行结果 #include #include static double T[200][200]; double fun_x( double x ) //被积函数 { return sqrt(x);} double Romberg( double a, double b, double Eps ) { int k=0; //用来记录把区间[a,b]2等分的次数 T[0][0] = (b-a)/2*( fun_x(a)+fun_x(b) ); do{ k++; double temp=0; for( int i=1; i<=pow(2,k-1); i++ ) { temp += fun_x(a+(2*i-1)*(b-a)/pow(2,k));} T[0][k] = 0.5*(T[0][k-1]+((b-a)/pow(2,k-1))*temp); for( int m=1; m<=k; m++ ) { T[m][k-m] = ( pow(4,m)*T[m-1][k-m+1]-T[m-1][k-m] )/( pow(4,m) - 1 );} }while(fabs( T[k][0]-T[k-1][0] )>=Eps); return T[k][0];} int main() { double a,b,Eps; cout<<"请输入积分上下限a,b及精度Eps:"; cin>>a>>b>>Eps; cout<<"所求积分的近似值为:"<return 0; }
第 15 页 共32页
本文来源:https://www.wddqw.com/doc/4b42bef3ae51f01dc281e53a580216fc700a531f.html