计算方法大作业——龙贝格积分

3 龙贝格积分

3.1 算法原理及程序框图

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

()2221

[]41

n n n n I f T T T T ≈+

-=-.

(11) 对于复化辛普森求积公式和复化科茨求积公式而言，也有类似的关系，如公式(12)和公式(13)。

()22221

[]41n n n n I f S S S S ≈+

-=- (12)

()22231

[]41

n n n n I f C C C C ≈+-=- (13)

通过对公式(11)~(13)做进一步分析，可得到公式(14)和公式(15)。

()221

41n n n n S T T T =+--

(14)

()2221

41

n n n n C S S S =+-- (15)

根据公式(14)和公式(15)表现出来的规律，令龙贝格积分为

()2231

41

n n n n R C C C =+-- (16)

其截断误差为c R h 8f (8)(η)，已经具有很高的精度。 龙贝格积分法是将区间[a , b ]逐次分半进行计算，因此，对已知函数f (x )在区间[a , b ]上的龙贝格积分法的计算公式的算法如下，程序框图如图13所示。

(1) 计算T 1：[]1()()2

b a

T f a f b -=+；

(2) 逐次计算T 2k +1：()1

211221

121,0,1,2222k

k k k k i b a b a T T f a i k +++=--⎛

⎫=++-= ⎪⎝⎭

∑；

(3) 逐次计算S 2k 、C 2k 和R 2k ：()()()111111222

222222322

22141141141k

k k k k k k k k k k k S T T T C S S S R C C C ++++++⎧

=+-⎪-⎪

⎪

=+-⎨-⎪

⎪

=+-⎪-⎩

；

(4) 若122k k R R ε+-<，则取[]12k I f R +≈；否则，继续计算，直到满足精度为止。

图13 龙贝格积分法程序框图

3.2 程序使用说明

龙贝格积分法的matlab 程序见附录4，该程序可以求解连续函数在特定区间上的近似积分值。运行程序按照命令窗口的提示输入相关变量直至得到结果。

开始运行程序时，命令窗口会出现提示“请输入被积函数表达式：f(x) =”，此时需

要输入被积函数的合法表达式，如被积函数为()sin x

f x x

=、()x f x e =+和

()2ln log f x x x =+则需要分别输入sin(x)/x 、exp(x)+sqrt(x)和log(x)+log2(x)；输入完成后按Enter 键进行下一步则出现提示“请输入积分区间下限：”，此时需要输入积分区间[a , b ]中的a 值，输入完成后继续按Enter 出现提示“请输入积分区间上限：”，根据提示输入b 的值，继续按Enter 出现提示“请输入积分误差：”，此时需要输入积分要求的精度，如要求的精度为10-12时输入1e-12即可。至此所有数据的输入已经完成，接着按Enter 键，程序会按照图13的顺序运行，直至给出最后的结果。

3.3 算例及计算结果

(1) 《数值分析》课本第186页的例题6.1.2，用龙贝格积分法计算积分

[]10sin x

I f dx x

=⎰

使误差不超过10-12。

使用龙贝格积分法的matlab 程序对该积分进行计算，运行结果如图14所示，最终可得积分的近似解为I [f ] ≈ 0.946083070367，其结果与课本的表6.1的结果吻合。

图14 龙贝格积分法计算例6.1.2的运行结果

(2) 《数值分析》课本第211页的计算实习题6.2，用龙贝格积分法计算下列积分，使计算结果尽可能准确。

③ 101

1dx x +⎰；

② 1

2

0ln(1x)

1dx x ++⎰；

③ 1

0ln(1)

x dx x

+⎰； ④

2

sin x

dx x

π

⎰

. 使用龙贝格积分法的matlab 程序对该积分进行计算，运行结果如图15~图17所示。

图15 积分①运行结果图16 积分②运行结果

图17积分③运行结果图18 积分④运行结果

由图15~图18运行结果可以看出，(2)中的四个积分结果最终分别为I[f]1 ≈0.822467033424、I[f]2 ≈ 0.272198261288、I[f]3 ≈ 0.693147180560和I[f]4 ≈ 1.370762168154。误差均为10-12。