龙贝格积分

来源：互联网发布：mac怎么下载ps破解版编辑：程序博客网时间：2024/05/13 10:53

算法原理和程序框图

龙贝格积分法是在复化梯形求积公式、复化辛普森求积公式和复化科茨求积公式关系的基础上，构造出的一种精度更高的数值积分方法。龙贝格积分为：

（1）

其截断误差为，已具有很高精度。龙贝格积分法是将区间[a, b]逐次分半进行计算。计算顺序为：计算T1，逐次计算T₂^k+1逐次计算S₂^k,C₂^k和R₂^k直到满足精度为止。因此，对已知函数f(x)在区间[a, b]上的龙贝格积分法的计算程序如下：

3.2程序使用说明

被积函数f（x），积分区间[a,b]和误差界都可根据需要在程序中直接修改，最大递推次数也可自己设定。给定好初值后直接运行就好。

3.3计算实例结果

选取第 211 页计算实习 6.2为例。通过设置断点逐个求解。被积函数为f(x)=1/(1+x)时，递推到5次后满足要求精度，所求积分I=0.6931472；被积函数为f(x)=ln(1+x)/(1+x^2)时，递推到4次后满足要求精度，所求积分I=0.2721983；被积函数为f(x)=ln(1+x)/x时，递推到5次后满足要求精度，所求积分I=0.8224670；被积函数为f(x)=sin(x)/x时，递推到4次后满足要求精度，所求积分I=0.9460831。

% ** 文件名:Romberg.m% % ** 创建人:**% % ** 日 期:2016.12.14% % ** 描 述:Romberg积分法计算定积分% % ** 函数：第 211 页 计算实习 6.2     % % ** 参考教材：《数值分析》李乃成，梅立泉，科学出版社%%%clear;clc% %被积函数为f(x)=1/(1+x);积分区间为[0,1];  保留7位有效数字%误差界eps=10^(-8)b=1;a=0;h=b-a;eps=10^(-8);kmax=30;%最大递推次数T1=h/2*(1/(a+1)+1/(b+1));S1=0;C1=0;C2=0;R1=0;R2=0;for k=1:kmax    h=(b-a)/2^k;    i=1:2^(k-1);    x=a+(2*i-1)*h;    fx=sum(1./(1+x));      T2=T1/2+fx*h ;    S2=T2+(T2-T1)/3;    if(k<3)        if k==2            C2=S2+(S2-S1)/15;        end    else        C2=S2+(S2-S1)/15;        R2=C2+(C2-C1)/63;        if abs(R2-R1)<eps            R2            R1            break;        end        R1=R2;    end    T1=T2;S1=S2;C1=C2;endfprintf('递推到%d次后满足要求精度\n',k);fprintf('所求积分I=%8.7f\n',R2);% 递推到5次后满足要求精度% 所求积分I=0.6931472%%clear;clc% %被积函数为f(x)=ln(1+x)/(1+x^2);积分区间为[0,1];  保留7位有效数字%误差界eps=10^(-8)b=1;a=0;h=b-a;eps=10^(-8);kmax=10;%最大递推次数T1=h/2*(log(1+a)/(1+a^2)+log(1+b)/(1+b^2));S1=0;C1=0;C2=0;R1=0;R2=0;for k=1:kmax    h=(b-a)/2^k;    i=1:2^(k-1);    x=a+(2*i-1)*h;    fx=sum(log(1+x)./(1+x.^2));      T2=T1/2+fx*h ;    S2=T2+(T2-T1)/3;    if(k<3)        if k==2            C2=S2+(S2-S1)/15;        end    else        C2=S2+(S2-S1)/15;        R2=C2+(C2-C1)/63;        if abs(R2-R1)<eps            R2            R1            break;        end        R1=R2;    end    T1=T2;S1=S2;C1=C2;endfprintf('递推到%d次后满足要求精度\n',k);fprintf('所求积分I=%8.7f\n',R2);% 递推到4次后满足要求精度% 所求积分I=0.2721983%%clear;clc% %被积函数为f(x)=ln(1+x)/x;积分区间为[0,1];  保留7位有效数字%误差界eps=10^(-8)b=1;a=0;h=b-a;eps=10^(-8);kmax=10;%最大递推次数T1=h/2*(1+(log(1+b))/b);%因为a不能取到0，在这里吧ln(1+x)/x在0处的极限值1放到里面S1=0;C1=0;C2=0;R1=0;R2=0;for k=1:kmax    h=(b-a)/2^k;    i=1:2^(k-1);    x=a+(2*i-1)*h;    fx=sum(log(1+x)./x);      T2=T1/2+fx*h ;    S2=T2+(T2-T1)/3;    if(k<3)        if k==2            C2=S2+(S2-S1)/15;        end    else        C2=S2+(S2-S1)/15;        R2=C2+(C2-C1)/63;        if abs(R2-R1)<eps            R2            R1            break;        end        R1=R2;    end    T1=T2;S1=S2;C1=C2;endfprintf('递推到%d次后满足要求精度\n',k);fprintf('所求积分I=%8.7f\n',R2);% % 递推到5次后满足要求精度% 所求积分I=0.8224670%%clear;clc% %被积函数为f(x)=sin(x)/x;积分区间为[0,pi/2];  保留7位有效数字%误差界eps=10^(-8)b=1;a=0;h=b-a;eps=10^(-8);kmax=10;%最大递推次数T1=h/2*(1+sin(b)/b);%因为a不能取到0，在这里吧sin(x)/x在0处的极限值1放到里面S1=0;C1=0;C2=0;R1=0;R2=0;for k=1:kmax    h=(b-a)/2^k;    i=1:2^(k-1);    x=a+(2*i-1)*h;    fx=sum(sin(x)./x);      T2=T1/2+fx*h ;    S2=T2+(T2-T1)/3;    if(k<3)        if k==2            C2=S2+(S2-S1)/15;        end    else        C2=S2+(S2-S1)/15;        R2=C2+(C2-C1)/63;        if abs(R2-R1)<eps            R2            R1            break;        end        R1=R2;    end    T1=T2;S1=S2;C1=C2;endfprintf('递推到%d次后满足要求精度\n',k);fprintf('所求积分I=%8.7f\n',R2);% 递推到4次后满足要求精度% 所求积分I=0.9460831

1 0