龙贝格积分的程序实现.docx

龙贝格积分的程序实现.docx

《龙贝格积分的程序实现.docx》由会员分享，可在线阅读，更多相关《龙贝格积分的程序实现.docx（3页珍藏版）》请在冰点文库上搜索。

计算方法实验报告3

【课题名称】

龙贝格积分的程序实现

【目的和意义】

函数变化有急有缓，为了照顾变化剧烈部分的误差，我们需要加密格点。

对于变化缓慢的部分，加密格点会造成计算的浪费。

以此我们介绍一种算法，可以自动在变化剧烈的地方加密格点计算，而变化缓慢的地方，则取稀疏的格点。

实际计算中，由于要事先给出一个合适的步长往往很困难，所以我们往往采用变步长的计算方案，即在步长逐步分半的过程中，反复利用复化求积公式进行计算，直到所求得的积分值满足精度要求为止。

在步长逐步分半过程中将粗糙的积分值逐步加工为精度较高的积分值，或者说将收敛缓慢的梯形值序列加工成收敛迅速的积分值序列。

这种加速方法称为龙贝格算法。

【计算公式】

设表示复化梯形求得的积分值，其下标是等分数，由此则有递推公式其中

，其中

由复化梯形公式的截断误差公式可得

，。

由此可知，。

这样导出的加速公式是辛普森公式：

同理可得。

由此便可得加速的算法：

龙贝格算法。

【龙贝格求积算法流程图】

定义被积函数f，积分上下限a，b和精度c

定义一个15×4的零矩阵，用于存放t值

先按照公式算出t（1,1）至t（1,3），以便循环

按公式计算其余t值

abs（t（k,4）-t（k-1,4））6

否

是

否

是

输出“不收敛”

输出近似值

k>=15

【龙贝格求积算法Matlab主程序】

function[t]=rbg（f,a,b,c）%定义龙贝格积分函数，f为待积函数，a与b为积分上下限，c为精度控制；

t=zeros（15,4）;%生成一零矩阵，用于存放t值；

t（1,1）=（b-a）/2*（f（a）+f（b））;%由于矩阵行列值均从1开始，所以将原本的t（0,0）记为t（1,1）,行列均加1；

fork=2:

4%先算出第一列的4个（包括t（1,1））值，以便程后面可以直接循环计算；

sum=0;

fori=1:

2^（k-2）

sum=sum+f（a+（2*i-1）*（b-a）/2^（k-1））;

end

t（k,1）=0.5*t（k-1,1）+（b-a）/2^（k-1）*sum;

fori=2:

k

t（k,i）=（4^（i-1）*t（k,i-1）-t（k-1,i-1））/（4^（i-1）-1）;

end

end

fork=5:

15%循环按照公式计算出t值；

sum=0;

fori=1:

2^（k-2）

sum=sum+f（a+（2*i-1）*（b-a）/2^（k-1））;

end

t（k,1）=0.5*t（k-1,1）+（b-a）/2^（k-1）*sum;

fori=2:

4

t（k,i）=（4^（i-1）*t（k,i-1）-t（k-1,i-1））/（4^（i-1）-1）;

end

ifk>6%可知最小二分次数，防止假收敛；

ifabs（t（k,4）-t（k-1,4））

disp（['答案',num2str（t（k,4））]）;

break;

end

end

end

ifk>=15

disp（['不收敛']）;%二分次数达15次仍不收敛；

end

【调用函数解题】

01sinxxdx

由于被积函数sin（x）/x在积分下限0时函数值难定，故取积分下限为10^（-60）。

调用函数解题如下：