Page 238 - Applied Numerical Methods Using MATLAB

TRAPEZOIDAL METHOD AND SIMPSON METHOD 227
N−1

f(a) + f(b)
I T 2 (a,b,h) = h + f(x k ) (5.6.1)
2
k=1
2
whose error is proportional to h as N times the error for one segment [Eq. (5.5.10)],
that is,
3
2
3
NO(h ) = (b − a)/h × O(h ) = O(h )
On the other hand, we have the numerical integration of f(x) over [a, b]by
Simpson’s rule (5.5.4) with an even number of segments N as
b N/2−1 x 2m+2

f(x) dx = f(x) dx
a x 2m
m=0
h
∼
= {(f 0 + 4f 1 + f 2 ) + (f 2 + 4f 3 + f 4 ) + ··· + (f N−2 + 4f N−1 + f N )}
3
N/2−1 N/2−1

h
I S4 (a,b,h) = f(a) + f(b) + 4 f(x 2m+1 ) + 2 f(x 2m ) (5.6.2)
3
m=0 m=1
N/2−1 N−1

h
= f(a) + f(b) + 2 f(x 2m+1 ) + f(x k )
3
m=0 k=1
4
whose error is proportional to h as N times the error for one segment [Eq. (5.5.11)],
that is,
4
5
5
(N/2)O(h ) = (b − a)/2h × O(h ) = O(h )
These two integration formulas by the trapezoidal rule and Simpson’s rule are
cast into the MATLAB routines “trpzds()”and “smpsns()”, respectively.
function INTf = trpzds(f,a,b,N)
%integral of f(x) over [a,b] by trapezoidal rule with N segments
if abs(b - a) < eps | N <= 0, INTf = 0; return; end
h = (b - a)/N; x = a +[0:N]*h; fx = feval(f,x); values of f for all nodes
INTf = h*((fx(1) + fx(N + 1))/2 + sum(fx(2:N))); %Eq.(5.6.1)
function INTf = smpsns(f,a,b,N,varargin)
%integral of f(x) over [a,b] by Simpson’s rule with N segments
if nargin < 4, N = 100; end
if abs(b - a)<1e-12 | N <= 0, INTf = 0; return; end
if mod(N,2) ~= 0, N=N+1;end %make N even
h = (b - a)/N; x = a + [0:N]*h; %the boundary nodes for N segments
fx = fevel(f,x,varargin{:}); %values of f for all nodes
fx(find(fx == inf)) = realmax; fx(find(fx == -inf)) = -realmax;
kodd = 2:2:N; keven = 3:2:N - 1; %the set of odd/even indices
INTf = h/3*(fx(1) + fx(N + 1)+4*sum(fx(kodd)) + 2*sum(fx(keven)));%Eq.(5.6.2)