利用Matlab编写龙贝格算法（romberg）求函数积分

 这次是我初次接触matlab，源于数学老师布置的一个作业：用龙贝格算法来计算函数的积分。

 

具体的计算原理，由于是数学的东西，不好打印，就不写了。主要把自己的代码贴下来慢慢理解。

 

一共写了两个文件。一个是romberg.m主要是写利用龙贝格算法，第二个是compute.m是调用之前写的接口

 

代码如下：

romberg.m

function [R,k,T]=romberg(fun,a,b,tol)
% 龙贝格(Romberg数值求解公式)
% author:
%   -gongwanlu
% inputs:
%   -fun：积分函数句柄
%   -a/b：积分上下限
%   -tol：积分误差
% Outputs:
%   -R：7阶精度Romberg积分值
%   -k：迭代次数
%   -T：整个迭代过程
%
% Example
% fun=@(x)4./(1+x^2);
% [R,k,T]=romberg(fun,0,1,1e-6)
%
k=0; % 迭代次数
n=1; % 区间划分个数
h=b-a;
T=h/2*(fun(a)+fun(b));
err=1;
while err>=tol
    k=k+1;
    h=h/2;
    tmp=0;
    for i=1:n
        tmp=tmp+fun(a+(2*i-1)*h);
    end
    T(k+1,1)=T(k)/2+h*tmp;
    for j=1:k
        T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);
    end
    n=n*2;
    err=abs(T(k+1,k+1)-T(k,k));
end
R=T(k+1,4);

 

compute.m

%计算（4/1+X^2）在0到1上面的积分
a = 0
b = 1
epsilon = 5e-6
f = @(x)4 ./ (1 + x^2);
y = romberg(f,a,b,epsilon) 


%后面是画出函数图像，注，不是积分函数图像。是被积函数图像
x = 0:0.01:1;
z = 4 ./ ( 1 + x.^2 );
plot(x,z),xlabel('x'),ylabel('y'),title(['Result = ',num2str(y)])

 

 画出的图像为，积分结果写在了图像上面。