数值666 (2)

【#文档大全网# 导语】以下是®文档大全网的小编为您整理的《数值666 (2)》，欢迎阅读！

龙贝格求积分算法

4.龙贝格求积分算法

4.1算法说明



生成JK的逼近表R(J,K)，并以R(J1,J1)为最终解来逼近积分













b

a

f(x)dxR(J,J)

逼近R(J,K)存在于一个特别的下三角矩阵中，第0列元素R(J,0)用基于2J个[a,b]子区间的连续梯形方法计算，然后利用龙贝格公式计算R(J,K)。当1KJ时，第J行的元素为













R(J,K)R(J,K1)

R(J,K1)R(J1,K1)



4K1

当|R(J,J)R(J1,J1)|tol时，程序在第(J1)行结束。

4.2龙贝格求积分算法流程图

图4-1 流程图



第 14 页 共32页


数值计算课程设计

4.3龙贝格求积分算法程序调试 求解积分方程果如图4-2：

9

1

xdx为例，对龙贝格求积分算法程序进行编译和链接，执行后得到结



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.wddqxz.cn/4b42bef3ae51f01dc281e53a580216fc700a531f.html