composite.m

来自「Integraton routines in matlab」· M 代码 · 共 62 行

%
% composite.m 
%
% Compute the weights of the Newton-Cotes quadrature quadrature
% on [0,1] with equidistant nodes.
%
%

% errror function
f = inline('exp(-x.^2)*2/sqrt(pi)', 'x');
a = 0; b = 2;
int = erf(b);

% bessel function
f = inline('cos( 8*sin(x) - x )/pi', 'x');
a = 0; b = pi;
int = bessel(1,8);

% 
f = inline('sin(x).^2.*cos(x).^4', 'x');
a = 0; b = 2*pi;
int = b/16;

% 
f = inline('4./(1+x.^2)', 'x');
a = 0; b = 1;
int = pi;

% composite trapezoidal rule
for n = 1:16
    % subinterval length
    h = (b-a)/n;
    % equidistant nodes
    x = (a:h:b);
    % function values at equidistant nodes
    fx = feval( f, x);
    % composite Trapezoidal rule
    T(n) = h*( 0.5*(fx(1)+fx(n+1)) + sum(fx(2:n)) );
end

% composite Simpson rule
for n = 1:18
    % subinterval length
    h = (b-a)/n;
    % equidistant nodes and midpoints
    x  = (a:h:b);
    % function values at equidistant nodes
    fx = feval( f, x);
    % composite Simpson rule
    S(n) = (h/6)*( fx(1)+fx(n+1) + 2*sum(fx(2:n)) );
    % midpoints
    xm = (a+h/2:h:b);
    % function values at midpoints
    fx = feval( f, xm);
    % composite Simpson rule
    S(n) = S(n) + (2*h/3)*sum(fx);
end

fprintf('    n  &       T(n)        &        S(n)       &       int \n' )
for n = 1:16
  fprintf('  %3d  &  %14.13f  &  %14.13f  & %14.13f \n', n+1, T(n), S(n), int )
end