注册 登录
编程论坛 Matlab

请教高手解决下面这题

ft4029928 发布于 2011-12-19 23:17, 572 次点击
用(梯形)龙贝格求积公式近似计算I=(积分上限1,下限0)2e(-x)/根号(x)dx 注:由于这里输不了公式,我用文字代替了,使其精度达到1e-5  
我运行时出错,找不到原因,请各位高手帮助我,非常感激!

程序代码:
function [val,M]=Romberg(f,ab,dalt)
% Romberg算法计 算积分
%  [val,M]=Romberg(f,ab,dalt)
%       val    返回积分值
%       M  T数表
%       f      被积函数,函数文件名
%       ab     积分区间  
%       dalt    精度
if nargin<3
   dalt=1e-7;                         %设置默认精度1e-7  
end
i=1; h=ab(2)-ab(1); t=0; j=0;
T(1,1)=(h/2)*(feval(f,ab)*ones(2,1));
while i<50                            %最大跌?0次
    i=i+1; h=h/2;
    T=[T zeros(i-1,1);zeros(1,i)];
    x=ab(1)+h:2*h:ab(2)-h;
    y=feval(f,x)*ones(size(x')); % feval(f,x)求f在x点值
    T(i,1)=T(i-1,1)/2+h*y;       % 细化区间,求得梯形值
    for t=2:i
        j=t-1;
        T(i,t)=(4^j*T(i,j)-T(i-1,j))/(4^j-1); %外推加速
    end
    if abs(T(i,i)-T(i-1,i-1))<=dalt     %控制精度
        break
    end
end
if nargout==2
   M=T;                                    %按需求返回T数表
end
val=T(i,i);                                %积分值

脚本文件:
y=inline('2*exp(-x).*(x.^(-1/2))');
[I,M]=Romberg(y,[0 1],1.0e-5)
0 回复
1