龙贝格积分算法实验

实验题目2 Romberg积分法

摘要

考虑积分

I(f)??欲求其近似值，可以采用如下公式： （复化）梯形公式 T?n?1i?0baf(x)dx

?2[f(x)?f(xihi?1)]

b?a2hf??(?) ??[a,b] 12n?1h （复化）辛卜生公式 S??[f(xi)?4f(x1)?f(xi?1)]

i?i?062 E??b?a?h?(4) E????f(?) ??[a,b ]180?2?n?1h （复化）柯特斯公式 C??[7f(xi)?32f(x1)?12f(x1)?

i?i?i?09042 32f(xi?344)?7f(xi?1)]

62(b?a)?h?(6) E????f(?) ??[a,b ]945?4?这里，梯形公式显得算法简单，具有如下递推关系

1hn?1T2n?Tn??f(x1)

i?22i?02因此，很容易实现从低阶的计算结果推算出高阶的近似值，而只需要花费较少的附加函数计算。但是，由于梯形公式收敛阶较低，收敛速度缓慢。所以，如何提高收敛速度，自然是人们极为关心的课题。为此，记T0,k为将区间[a,b]进行2等份的复化梯形积分结果，T1,k为将区间[a,b]进行2等份的复化辛卜生积分结果，T2,k为将区间[a,b]进行2等份的复化柯特斯积分结果。根据李查逊（Richardson）外推加速方法，可得到

kkkTm,k?4mTm?1,k?1?Tm?1,k4m?1?k?0,1,2,???? m?0,1,2,??? 可以证明，如果f(x)充分光滑，则有

limTm,k?I(f) （m固定）

k?? limTm,0?I(f)

m??这是一个收敛速度更快的一个数值求积公式，我们称为龙贝格积分法。 该方法的计算可按下表进行

T0,0 T0, 1 T0, 2 … T0,m T1,0 T1, 1,m?1 1 … T T2,0 … T2,m?2 … … Tm,0

很明显，龙贝格计算过程在计算机上实现时，只需开辟一个一维数组，即每次计算的结果

Tm,k，可存放在T0,k位置上，其最终结果Tm,0是存放在T0,0位置上。

前言

利用龙贝格积分法计算积分I(f)??baf(x)dx

程序设计流程：

1．准备初值，计算

T0,0?且k?0（k为等份次数）

2．按梯形公式的递推关系，计算

a?b[f(a)?f(b)] 21b?a2?1b?a1T0,k?1?T0,k?k?1?f(a?k(i?))

2222i?0 3．按龙贝格公式计算加速值

kT0,k?m?Tm,k?m?4mTm?1,k?1?m?Tm?1,k?m4?1mm?0,1,2,?,k

4．精度控制。对给定的精度?，若

|Tm,0?Tm?1,0|??

则终止计算，并取T0,s?Tm,s作为所求结果；否则k?k?1，重复2~4步，直到满足精度为止。

问题1

（1） 程序运行如下：

I = Romberginterg(inline('x.^2.*exp(x)'),0,1,25,1e-6) I = 0.7183

（2） 程序运行如下：

I = Romberginterg(inline('exp(x).*sin(x)'),1,3,25,1e-6) I = 10.9502

（3） 程序运行如下：

I = Romberginterg(inline('4./(1+x.^2)'),0,1,25,1e-6) I = 3.1416

（4） 程序运行如下：

I = Romberginterg(inline('1./(1+x)'),0,1,25,1e-6) I = 0.6931

问题2

（1） 程序运行如下：

I = Romberginterg(inline('cos(x)./sqrt(1-x)'),0,1,25,1e-6) I = NaN

因为函数在x=0处出现了0/0的情况，极限为1，所以Matlab的结果显示为NaN非数，解决方法是把下限0改为一个小数eps。

I = Romberginterg(inline('sin(x)./x'),eps,1,25,1e-6) I = 0.9461 （2） 程序运行如下：

I = Romberginterg(inline('cos(x)./sqrt(1-x)'),0,1,25,1e-6) I = NaN

与（1）的原因相同，函数在x=1处出现了0/0的情况，结果为无穷，此时应选择变换

，

将积分变为，再进行变换，将积分变为

，变换后输入命令：

I = RombergInterg(inline('2*sin(x).*cos((sin(x)).^2)'),0,pi/2,25,1e-6) I = 1.4996

（3） 程序运行如下：

I = Romberginterg(inline('cos(x)./sqrt(x)'),0,1,25,1e-6)

I = NaN

函数在x=0处出现了1/0的情况，结果为无穷。先将积分变为

，再做

变换，

I = 2*Romberginterg(inline('cos(x.^2)'),0,1,25,1e-6) I = 1.8090

（4） 程序运行如下：

I = Romberginterg(inline('x.*sin(x)./(1-x.^2)'),-1,1,25,1e-6) I = NaN

函数在x=1,-1处出现了sin1/0的情况，结果为无穷。被积函数为偶函数，做变换

，

积分变为

I = 2*Romberginterg(inline('sin(x).*sin(sin(x))'),0,pi/2,25,1e-6) I = 1.3825

或者使用Gauss-Chebyshev求积公式：

I = GaussInterg(inline('x.*sin(x)'),'Chebyshev',-1,1,1e-6) I = 1.3825

由于Gauss-Chebyshev求积公式求的是。

在

的积分因此此处的

为

问题3

（1） 程序运行如下：

I = Romberginterg(inline('exp(-x.^2)'),-10,10,25,1e-6) I = 1.7725

由于积分区间无限，而函数在[-10,10]外很小，可以用函数在[-10,10]上的积分近似这个无穷积分。或者使用Gauss-Hermite求积公式：

I = GaussInterg(inline('x+1-x'),'Hermite',0,0,1e-6) I = 1.7725

由于Gauss-Hermite求积公式求的是（2） 程序运行如下：

先做变换

，则原积分变为

在

的积分因此此处的

为1。

I = 2*Romberginterg(inline('1./(1+x.^2)'),0,1,25,1e-6)

I = 1.5708

（3） 程序运行如下：

使用Gauss-Hermite求积公式：

I = GaussInterg(inline('cos(x).^3'),'Hermite',0,0,1e-6) I = 1.0820

由于Gauss-Hermite求积公式求的是

在

的积分因此此处的

为

。

（4） 程序运行如下：

使用Gauss-Laguerre求积公式：

I = GaussInterg(inline('sin(x).^2'),'Laguerre',0,0,1e-6) I = 0.4000

由于Gauss-Laguerre求积公式求的是

在

的积分因此此处的

为

。

使用的函数

function I = RombergInterg(fun, a, b, npanel, tol, flag) % RombergInterg 用Romberg方法求积分 %

% Synopsis: I = RombergInteg(fun,a,b,n,tol) %

% Input: fun = (string) 被积函数的函数名 % a, b = 积分下限和积分上限

% npanel = (optional) 将积分区间平分的段数，默认为25 % tol = (optional) 计算误差上限，默认为5e-9

% flag = (optional) 是否显示Romberg-T数表，默认为0为不显示 %

% Output: I = 通过Romberg方法求积分的近似值 if nargin < 4

npanel = 25; end

if nargin < 5 tol = 5e-9; end

if nargin < 6 flag = 0; end

T(1,1) = TrapezoidInteg(fun, a, b, npanel); %T0(h) = T(h)

err = 1; %初始化误差值

m = 2; %初始化行计算的行序号

while err >= tol

T(m,1) = TrapezoidInteg(fun, a, b, 2^m*npanel); %计算第m行T0(h/2^m) T(m,m) = 0; %将矩阵T变成m*m

for n = 2:m

T(m,n) = ( 4^(n-1)*T(m,n-1) - T(m-1,n-1)) / (4^(n-1) - 1); %Tm(h/2^k)与Tm-1(h/2^(k+1))和Tm-1(h/2^k)的递推关系 end

err = abs( T(m,m) - T(m-1,m-1) ); %计算误差值

m = m + 1; %计算下一行 end

I = T(m-1,m-1);

if flag ~= 0 disp(T); end

function I = TrapezoidInteg(fun, a, b, npanel) % TrapezoidInteg 用复化梯形公式求积分 %

% Synopsis: I = TrapezoidInteg(fun,a,b,n) %

% Input: fun = (string) 被积函数的函数名 % a, b = 积分下限和积分上限

% npanel = (optional) 将积分区间平分的段数，默认为25 %

% Output: I = 通过复化梯形公式求积分的近似值

if nargin < 4

npanel = 25; end

nnode = npanel + 1; %节点数 = 段数 + 1 h = (b-a)/(nnode-1); %步长

x = a:h:b; %将积分区间分段 f = feval(fun,x); %求节点处被积函数的值

I = h * ( 0.5*f(1) + sum(f(2:nnode-1)) + 0.5*f(nnode) );

function I = GaussInterg(fun, type, a, b, tol)

% GaussInterg 用Gauss型求积公式求积分,具体形式由使用者选取 %

% Synopsis: I = GaussInterg(fun, type, a, b)

% I = GaussInterg(fun, type, a, b, tol) %

% Input: fun = (string) 被积函数的函数名

% type = (string) 具体Gauss求积公式的形式

% a,b = 积分上下限，Laguerre只计算0到inf，Hermite只计算-inf到inf，所以a,b对这两种形式无效

% tol = (optional) 误差容忍限度，默认为5e-5 %

% Output: I = 通过Gauss型求积公式求积分的近似值 if nargin < 5 tol = 5e-5; end

n = 7; %默认从7节点多项式开始计算 IOld = 1; %初始化IOld为1 err = 1; %初始化误差为1 while err >= tol switch type

case 'Legendre' %计算无奇点被积函数在-1到1的积分 I = GaussLegendreInterg(fun, a, b, n);

case 'Chebyshev' %计算被积函数形如 f(x)/sqrt(1-x^2)在-1到1的积分 I = GaussChebyshevInterg(fun, a, b, n);

case 'Laguerre' %计算被积函数形如 exp(-x)*f(x)的在0到inf的积分 I = GaussLaguerreInterg(fun, n);

case 'Hermite' %计算被积函数形如 exp(-x^2)*f(x)的在-inf到inf的积分 I = GaussHermiteInterg(fun, n); otherwise

error('No such type!'); end

err = abs(I - IOld); %计算误差

IOld = I; %把IOld赋值为I进行下次迭代 n = n+1; %多项式节点递增 end

function I = GaussLegendreInterg(fun, a, b, n)

% GaussLegendreInterg 用Gauss-Legendre型求积公式计算无奇点被积函数在-1到1的积分 %

% Synopsis: I = GaussLegendreInterg(fun, a, b) %

% Input: fun = (string) 被积函数的函数名 % a,b = 积分下限和积分上限

% n = (optional) 高斯节点的个数，默认为7 %

% Output: I = 通过Gauss型求积公式求积分的近似值 if nargin < 4 n = 7; end

if round(n) ~= n | n < 1

error('Gauss节点个数必须为正整数'); end

p = LegendreIter(n); %构造n阶勒让德多项式 p = p/polyval(p,1); %使Pn(1) = 1;

for i = 1:n+1 %构造p的导数p1 p1(i) = (n-i+1)*p(i); end

p1(n+1) = [];

r = roots(p); %计算勒让德多项式的根

L = polyval(p1,r); %计算与勒让德多项式的根匹配的高斯积分系数 A = 2./(1-r.^2)./L.^2;

xi = ( (b-a)*r + (a+b) ) / 2; %找出勒让德多项式的根在[a,b]上对应的数值xi yi = feval(fun, xi); %计算f(xi)

I = (b-a)/2 * (A'*yi); %I = sum(Ai*yi),前面的系数是积分区间改变是增加的常数

function I = GaussChebyshevInterg(fun, a, b, n)

% GaussChebyshevInterg 用Gauss-Chebyshev型求积公式求积分计算被积函数形如 f(x)/sqrt(1-x^2)在-1到1的积分 %

% Synopsis: I = GaussChebyshevInterg(fun, a, b) %

% Input: fun = (string) 被积函数的函数名 % a,b = 积分下限和积分上限

% n = (optional) 高斯节点的个数，默认为7 %

% Output: I = 通过Gauss型求积公式求积分的近似值 if nargin < 4 n = 7; end

if round(n) ~= n | n < 1

error('Gauss节点个数必须为正整数'); end

p = ChebyshevIter(n); %构造n阶切比雪夫多项式

A = pi/n * ones(n,1); %计算与切比雪夫多项式的根匹配的高斯积分系数

r = roots(p); %计算切比雪夫多项式的根

xi = ( (b-a)*r + (a+b) ) / 2; %找出切比雪夫多项式的根在[a,b]上对应的数值xi yi = feval(fun, xi); %计算f(xi)

I = (b-a)/2 * (A'*yi); %I = sum(Ai*yi),前面的系数是积分区间改变是增加的常数

function I = GaussLaguerreInterg(fun, n)

% GaussLaguerreInterg 用Gauss-Laguerre型求积公式求积分计算被积函数形如 exp(-x)*f(x)的在0到inf的积分 %

% Synopsis: I = GaussLaguerreInterg(fun)

% I = GaussLaguerreInterg(fun, n) %

% Input: fun = (string) 被积函数的函数名

% n = (optional) 高斯节点的个数，默认为7 %

% Output: I = 通过Gauss型求积公式求积分的近似值 if nargin < 2 n = 7; end

if round(n) ~= n | n < 1

error('Gauss节点个数必须为正整数'); end

p = LaguerreIter(n); %构造n阶拉盖尔多项式 p1 = LaguerreIter(n+1); %构造n+1阶拉盖尔多项式

xi = roots(p); %计算拉盖尔多项式的根

L = polyval(p1,xi); %计算与拉盖尔多项式的根匹配的高斯积分系数 A = xi./((n+1)^2*L.^2);

yi = feval(fun, xi); %计算f(xi)

I = A'*yi; %I = sum(Ai*yi),前面的系数是积分区间改变是增加的常数

function I = GaussHermiteInterg(fun, n)

% GaussHermiteInterg 用Gauss-Hermite型求积公式求积分计算被积函数形如 exp(-x^2)*f(x)的在-inf到inf的积分 %

% Synopsis: I = GaussHermiteInterg(fun)

% I = GaussHermiteInterg(fun, n) %

% Input: fun = (string) 被积函数的函数名

% n = (optional) 高斯节点的个数，默认为7 %

% Output: I = 通过Gauss型求积公式求积分的近似值 if nargin < 2 n = 7; end

if round(n) ~= n | n < 1

error('Gauss节点个数必须为正整数'); end

p = HermiteIter(n); %构造n阶埃尔米特多项式 p1 = HermiteIter(n-1); %构造n-1阶埃尔米特多项式

xi = roots(p); %计算埃尔米特多项式的根

L = polyval(p1(1:n),xi); %计算与埃尔米特多项式的根匹配的高斯积分系数 A = 2^(n-1)*sqrt(pi)*factorial(n) ./ (n^2 * L.^2);

yi = feval(fun, xi); %计算f(xi)

I = A'*yi; %I = sum(Ai*yi),前面的系数是积分区间改变是增加的常数

function p = LegendreIter(n)

% LegendreIter 用递推的方法计算n次勒让德多项式的系数向量 Pn+2(x) = (2*i+3)/(i+2) * x*Pn+1(x) - (i+1)/(i+2) * Pn(x) %

% Synopsis: p = LegendreIter(n) %

% Input: n = 勒让德多项式的次数 %

% Output: p = n次勒让德多项式的系数向量 if round(n) ~= n | n < 0

error('n必须是一个非负整数'); end

if n == 0 %P0(x) = 1 p = 1; return;

elseif n == 1 %P1(x) = x p = [1 0]; return; end

pBk = 1; %初始化三项递推公式后项为P0 pMid = [1 0]; %初始化三项递推公式中项为P1 for i = 0:n-2

pMidCal = zeros(1,i+3); %构造用于计算的x*Pn+1 pMidCal(1:i+2) = pMid;

pBkCal = zeros(1,i+3); %构造用于计算的Pn pBkCal(3:i+3) = pBk;

pFwd = (2*i+3)/(i+2) * pMidCal - (i+1)/(i+2) * pBkCal; %勒让德多项式三项递推公式Pn+2(x) = (2*i+3)/(i+2) * x*Pn+1(x) - (i+1)/(i+2) * Pn(x)

pBk = pMid; %把中项变为后项进行下次迭代 pMid = pFwd; %把前项变为中项进行下次迭代 end

p = pFwd;

function p = ChebyshevIter(n)

% ChebyshevIter 用递推的方法计算n次勒让德-切比雪夫多项式的系数向量 Tn+2(x) = 2*x*Tn+1(x) - Tn(x) %

% Synopsis: p = ChebyshevIter(n) %

% Input: n = 勒让德-切比雪夫多项式的次数 %

% Output: p = n次勒让德-切比雪夫多项式的系数向量 if round(n) ~= n | n < 0

error('n必须是一个非负整数'); end

if n == 0 %T0(x) = 1 p = 1; return;

elseif n == 1 %T1(x) = x p = [1 0]; return; end

pBk = 1; %初始化三项递推公式后项为T0 pMid = [1 0]; %初始化三项递推公式中项为T1 for i = 0:n-2

pMidCal = zeros(1,i+3); %构造用于计算的x*Tn+1 pMidCal(1:i+2) = pMid;

pBkCal = zeros(1,i+3); %构造用于计算的Pn pBkCal(3:i+3) = pBk;

pFwd = 2*pMidCal - pBkCal; %勒让德-切比雪夫多项式三项递推公式Tn+2(x) = 2*x*Tn+1(x) - Tn(x)

pBk = pMid; %把中项变为后项进行下次迭代 pMid = pFwd; %把前项变为中项进行下次迭代 end

p = pFwd;

function p = LaguerreIter(n)

% LaguerreIter 用递推的方法计算n次拉盖尔多项式的系数向量 Ln+2(x) = (2*n+3-x)*Ln+1(x) - (n+1)*Ln(x) %

% Synopsis: p = LaguerreIter(n) %

% Input: n = 拉盖尔多项式的次数 %

% Output: p = n次拉盖尔多项式的系数向量 if round(n) ~= n | n < 0

error('n必须是一个非负整数'); end

if n == 0 %L0(x) = 1 p = 1; return;

elseif n == 1 %L1(x) = -x+1 p = [-1 1]; return; end

pBk = 1; %初始化三项递推公式后项为L0 pMid = [-1 1]; %初始化三项递推公式中项为L1 for i = 0:n-2

pMidCal1 = zeros(1,i+3); %构造用于计算的x*Ln+1(x) pMidCal1(1:i+2) = pMid;

pMidCal2 = zeros(1,i+3); %构造用于计算的Ln+1(x) pMidCal2(2:i+3) = pMid;

pBkCal = zeros(1,i+3); %构造用于计算的Ln(x) pBkCal(3:i+3) = pBk;

pFwd =( (2*i+3)*pMidCal2 - pMidCal1 - (i+1)*pBkCal )/ (i+2); %拉盖尔多项式三项递推公式Ln+2(x) = (2*n+3-x)*Ln+1(x) - (n+1)^2*Ln(x)

pBk = pMid; %把中项变为后项进行下次迭代 pMid = pFwd; %把前项变为中项进行下次迭代 end

p = pFwd;

function p = HermiteIter(n)

% HermiteIter 用递推的方法计算n次埃尔米特多项式的系数向量 Hn+2(x) = 2*x*Hn+1(x) - 2*(n+1)*Hn(x) %

% Synopsis: p = HermiteIter(n) %

% Input: n = 埃尔米特多项式的次数 %

% Output: p = n次埃尔米特多项式的系数向量 if round(n) ~= n | n < 0

error('n必须是一个非负整数'); end

if n == 0 %H0(x) = 1 p = 1; return;

elseif n == 1 %H1(x) = 2*x p = [2 0]; return; end

pBk = 1; %初始化三项递推公式后项为L0 pMid = [2 0]; %初始化三项递推公式中项为L1 for i = 0:n-2

pMidCal = zeros(1,i+3); %构造用于计算的x*Hn+1(x) pMidCal(1:i+2) = pMid;

pBkCal = zeros(1,i+3); %构造用于计算的Hn(x) pBkCal(3:i+3) = pBk;

pFwd =2*pMidCal - 2*(i+1)*pBkCal; %埃尔米特多项式三项递推公式Hn+2(x) = 2*x*Hn+1(x) - 2*(n+1)*Hn(x)

pBk = pMid; %把中项变为后项进行下次迭代 pMid = pFwd; %把前项变为中项进行下次迭代 end

p = pFwd;

思考题

（1） 输入的参数N越大，在有限区间上分段越小，由复化梯形公式的误差

，计算精度越高。

（2） 二分次数越多，计算精度越高。因为每一次二分，N变为原来的二倍，由复化梯形

公式的误差

，二分后误差约为二分前的四分之一。

（3） 有普遍性，被积函数在积分区间内存在奇点。若在奇点上极限存在，则可以利用小

数—机器Epsilon将奇点从积分区间去除，如问题2（1）。若奇点上极限不存在，则需要根据具体情况进行积分变换使被积函数在积分区间上不存在奇点。若求的是

在

的积分，则可以使用Gauss-Chebyshev求积公式。

在

的积分可用

（4） 有普遍性，积分区间均为无限区间。若求的是

Gauss-Hermite求积公式，若求的是

在

的积分可用Gauss-Laguerre求积

公式，其他形式根据具体情况进行积分变换使积分区间变为有限区间。

本文来源：https://www.bwwdw.com/article/v8e8.html