龍貝格演算法matlab

⑴ 用龍貝格法求積分

先用另外2種方法。

format long
%【1】精確值。符號積分
it=int('(2/sqrt(pi))*exp(-x)',0,1)
Accurate=eval(it)
y=inline('(2/sqrt(pi))*exp(-x)')
%【2】Simpson方法
Simpson=quad(y,0,1)
delta=Simpson-Accurate

結果：
Accurate = 0.713271669674918

y = Inline function:
y(x) = (2/sqrt(pi))*exp(-x)

Simpson = 0.713271671228492

delta = 1.553574158208448e-009

【3】從網上找到一個，存為romberg.m
%=================
function R = romberg(f, a, b, n)
format long
% ROMBERG -- Compute Romberg table integral approximation.
%
% SYNOPSIS:
% R = romberg(f, a, b, n)
%
% DESCRIPTION:
% Computes the complete Romberg table approximation to the integral
%
% / b
% I = | f(x) dx
% / a
%
% PARAMETERS:
% f - The integrand. Assumed to be a function callable as
% y = f(x)
% with `x' in [a, b].
% a - Left integration interval endpoint.
% b - Right integration interval endpoint.
% n - Maximum level in Romberg table.
%
% RETURNS:
% R - Romberg table. Represented as an (n+1)-by-(n+1) lower
% triangular matrix of integral approximations.
%
% SEE ALSO:
% TRAPZ, QUAD, QUADL.

% NOTE: all indices adjusted for MATLAB's one-based indexing scheme.

% Pre-allocate the Romberg table. Avoids subsequent re-allocation which
% is often very costly.
R = zeros([n + 1, n + 1]);

% Initial approximation. Single interval trapezoidal rule.
R(0+1, 0+1) = (b - a) / 2 * (feval(f, a) + feval(f, b));

% First column of Romberg table. Increasingly accurate trapezoidal
% approximations.
for i = 1 : n,
h = (b - a) / 2^i;

s = 0;
for k = 1 : 2^(i-1),
s = s + feval(f, a + (2*k - 1)*h);
end

R(i+1, 0+1) = R(i-1+1, 0+1)/2 + h*s;
end

% Richardson extrapolation gives remainder of Romberg table.
%
% Note: The table is computed by columns rather than the more
% traditional row version. The reason is that this prevents frequent
% (and needless) re-computation of the `fac' quantity.
%
% Moreover, MATLAB matrices internally use ``column major'' ordering so
% this version is less harmful to computer memory cache systems. This
% reason is an implementational detail, though, and less important in
% introctory courses such as MA2501.
for j = 1 : n,
fac = 1 / (4^j - 1);
for m = j : n,
R(m+1, j+1) = R(m+1, j-1+1) + fac*(R(m+1, j-1+1) - R(m-1+1, j-1+1));
end
end

function ff=f(x)
ff=2/sqrt(pi)*exp(-x);
%=================

運行：
>> R=romberg('f', 0, 1, 5)

R =

0.771743332258054 0 0 0 0 0
0.728069946441243 0.713512151168973 0 0 0 0
0.716982762290904 0.713287034240791 0.713272026445579 0 0 0
0.714200167058928 0.713272635314936 0.713271675386546 0.713271669814180 0 0
0.713503839348432 0.713271730111600 0.713271669764711 0.713271669675476 0.713271669674932 0
0.713329714927254 0.713271673453528 0.713271669676323 0.713271669674920 0.713271669674918 0.713271669674918

⑵ 用matlab，龍貝格演算法計算∫（0到1）[x/（4+x²）]dx的近似值。 求程序代碼！！！

%龍貝格求積演算法function I=romberg(a,b)h=b-a;T(1)=h/2*(fun(a)+fun(b));m=1;while 1 h=h/2; S(1)=1/2*T(1)+h*sumf(2^(m-1),a,h); for j=1:m S(j+1)=S(j)+(S(j)-T(j))/(4^j-1); end if abs(S(m+1)-T(m))<1e-6 break; end T=S;m=m+1;endI=S(m+1);end
function f=sumf(m,a,h)for j=1:m y(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
>>

⑶ 龍貝格演算法 橢圓周長 請幫忙修改一下程序~！！！！！！

1、在for(i=1;i<=pow(2,k-1);i++)之前加t=0;否則t重復累加；

2、for(i=1;i<pow(2,k-1);i++)改成for(i=1;i<=pow(2,k-1);i++)，i需執行到i=pow(2,k-1)；

3、調換for(k=0;k<20;k++)

{

for(m=1;m<4;m++)

的順序，改為

for(m=1;m<4;m++)

{

for(k=0;k<21-m;k++)否則在算T[k][m]時，T[k+1][m-1]還未知

4、並把k<20改成k<21-m；因為辛普森序列比梯形序列少一項，科茨序列比辛普森序列少一項，龍貝格序列比科茨序列少一項

話說哦。k取20是不太大了點啊，2的20次方可是很大的咯。其實k取6應該就能達到精度要求了吧。