用龙贝格法求积分被积函数为（2/sqrt(pi)）*exp（-x）区间为（0，1），要求误差不超过。（计算取6位小数）

Question

用龙贝格法求积分被积函数为（2/sqrt(pi)）*exp（-x）区间为（0，1） ，要求误差不超过 。（计算取6位小数）

Answer 1

先用另外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 % introductory 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