用MATLAB怎么编写Romberg方法

function [y,n,yiter]=RombergInt(fun,a,b,epsilon)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%Author:wacs5

%%Email :wacs5@126.com

%%Date :20100506

%%Function :Romberg积分

%%Argument :fun是被积函数；a,b是积分上下限；epsilon是误差容限

%%Return:y为积分结果；将[a,b]平均n等分；yiter迭代的积分值

%%Called Mtb:abs,feval,mfilename,nargin,sum,warning

%%Example :[y,n,yiter]=RombergInt('sin',0,pi/2,1e-10)

%%NOTICE:龙贝格积分返回的yiter含有计算过程的所有积分值

%%NOTICE:(个数为sum(1:log2(n)+1))

%%Reference :徐萃薇,孙绳武.计算方法引论(第二版).高教出版社 P117

%%See also :http://www.matlabsky.com/thread-7208-1-1.html

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (nargin<4)

epsilon=1e-6

end

if (nargin<3)

warning(['usage:',mfilename,'(fun,a,b[,epsilon])'])

end

n=1

h=(b-a)/2

y0=h*(feval(fun,a)+feval(fun,b))

yiter=y0

k=1

%%%%%%Variable%%%%%%%%%%%%%

%%y0:上次迭代的斜线

%%y1:这次迭代的斜线

%%yiter:保存所有的迭代值

while 1

y1(1)=y0(1)/2+h*sum(feval(fun,a+(1:2:2*n-1)*h))

mm=1

for m=1:k

mm=mm*4

y1(m+1)=(mm*y1(m)-y0(m))/(mm-1)

end

if abs(y1(end)-y0(end))<epsilon

break

end

n=n+n

h=h/2

k=k+1

y0=y1

yiter=[yiter,y1]

end

y=y1(end)

a=0

b=1

n=1

k=1

e=0.5*10^(-6)

T=zeros(10,10)

h=(b-a)/2

t=h*(1+sin(1))

T(1,1)=t

for j=1:9

F=0

for i=1:n

F=F+sin(a+(2*i-1)*h)/(a+(2*i-1)*h)

end

for i=1:k

if T(1,i+1)==0

T(1,i+1)=T(1,i)/2+h*F

end

end

for m=1:k

if T(m+1,k-m+1)==0

T(m+1,k-m+1)=(4^m*T(m,k-m+2)-T(m,k-m+1))/(4^m-1)

end

end

if abs(T(m+1,1)-T(m,1))<e

I=T(m+1,1)

break

else

h=h/2

n=2*n

k=k+1

end

end

最后的 I 就是函数值，K就是满足精度时的M值

%龙贝格求积算法function I=romberg(a,b)h=b-aT(1)=h/2*(fun(a)+fun(b))m=1while 1h=h/2 S(1)=1/2*T(1)+h*sumf(2^(m-1),a,h) for j=1:mS(j+1)=S(j)+(S(j)-T(j))/(4^j-1) endif abs(S(m+1)-T(m))<1e-6break endT=Sm=m+1endI=S(m+1)end

function f=sumf(m,a,h)for j=1:my(j)=fun(a+(2*j-1)*h)endf=sum(y)end

function f=fun(x)f=x/(4+x^2)end结果：

>>I=romberg(0,1)

I =

0.111571775646293

>>

---