龙贝格算法

第1章 龙贝格算法
研究问题阐述
计算二重积分
⎰⎰
=
b a
d
c
dxdy y x f I ),(------------------------------------(1-1)
对工科生要求计算如下3个二重积分：
算例1：11
10
()1I x y dxdy =
+=⎰⎰
算例2：/2
200
()I x y dxdy πππ=+=⎰⎰
-----------------------(1-2)
算例3： 1
1
30
sin()x y I dxdy x y
+=
+⎰⎰
通常将积分区间等分数依次取2k ，并得到相应的递推式子如（2-5）所示，相应对到过程不再累述详情可以参考教材。

111
2
22
1[()()],21[(21)]222k k k k k i b a T f a f b b a b a T T f a i --=-⎧=+⎪⎪
⎨--⎪=++-⎪⎩
∑----------------------------(2-5) 最终得到计算机计算步骤如下：
(1)计算初值T1； (2)K 1；
(3)计算新的T2，一直迭代计算2k
T ;
(4)精度控制：如果1
22
||k
k T T ε--<，则停止计算，并将2k T 作为积分的近似数值，否则
循环将不停止。

对第四点，通常可以加入一个循环最高次数，防止程序陷入死循环。

function f=fxy(x,y,w) if w==1
f=x+y; %函数1
elseif w==2
f=sin(abs(x-y)); %函数2
elseif w==3
f=sin(x+y)/(x+y); %函数3
if x+y==0
f=1;
end;
end;
function gy=gy(y,a,b,w,error)
k=1;lock=0;count=0;
T1=(b-a)*(fxy(a,y,w)+fxy(b,y,w))/2;
while lock==0&&count<100 %限定最多迭代20次
T2=T2kx(y,T1,a,b,k,w);
if abs(T2-T1)<error
gy=T2;
lock=1;
else
T1=T2;
count=count+1;
k=k+1;
end;
end
% count
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%
%%%%%%%%%%%%%%%%%%%% 龙贝格算法%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%% 作者：XXXX %%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%% 时间：2012/4/16 %%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%
clear all; clc;
%% 详情见报告Romberg公式介绍部分
clear all;clc;
which=1; %表示计算的函数f1,f2,f3
%积分区间
a=0;
if which==1||which==3
b1=1;b2=1; %函数1和2的区间else
b1=pi;b2=pi/2; %函数2的区间
end
%% 计算主程序
k=1;lock=0;count=0; %辅助量
error=0.0001; %计算精度
% gy=gy(y,a,b,w,error)
T1=(b2-a)*(gy(a,a,b1,which,error)+gy(b2,a,b1,which,error))/2; while lock==0&&count<100 %限定最多迭代20次T2=T2ky(T1,a,b1,b2,k,which,error);
if abs(T2-T1)<error
I=T2;
lock=1;
else
T1=T2;
function T2kx=T2kx(y,T2kx_1,a,b,k,w)
temp1=0;
for i=1:2^(k-1)
temp1=temp1+fxy(a+(2*i-1)*(b-a)/2^k,y,w);
end
T2kx=0.5*T2kx_1+temp1*(b-a)/2^k;
function T2ky=T2ky(T2ky_1,a,b1,b2,k,w,error)
temp1=0;
for i=1:2^(k-1)
temp1=temp1+gy(a+(2*i-1)*(b2-a)/2^k,a,b1,w,error);
end
T2ky=0.5*T2ky_1+temp1*(b2-a)/2^k;。