Simpson's Rule and Trapezoidal Rule of Numerical Integration - Matlab Scripts

Simpson's Rule and Trapezoidal Rule of Numerical Integration calculates the area under a function by breaking the function up into many smaller areas which are easier to calculate.

Trapezoidal Rule

>> trapezoidal('sin(x)/x',1,3,8)

ans =

   0.902337806742469

I = (((3-1)/8)/2)[f(1)+2{f(1.25)+f(1.5)+f(1.75)+f(2)+f(2.25)+f(2.5)+f(2.75)}+f(3)]

  = 0.125[0.841471+2(3.165097)+0.047040)

  = 0.902338

Trapezoidal Rule for numerical integration using trapazoids to calculate the area under the function. Each trapazoid models a linear line between two points on the function and calculates the area under the line.

Simpson's Rule

>> simprule('sin(x)/x',1,3,8)

ans =

   0.902568788567005

Instead of using linear lines to model the function, Simpson's Rule calculates the area by modelling a polynomial to the function. This polynomial is precise for other polynomials of degree 2 and 3.

I = (((3-1)/8)/3)[f(1)+4{f(1.25)+f(1.75)+f(2.25)+f(2.75)}+2{f(1.5)+f(2)+f(2.5)}+f(3)]

  = (0.25/3)[0.841471+4(1.806062)+2(1.359035)+0.047040)

  = 0.902569

Trapezoidal Rule Matlab Script

function I = trapezoidal(f_str, a, b, n)
%TRAPEZOIDAL Trapezoidal Rule integration.
% I = TRAPEZOIDAL(F_STR, A, B, N) returns the Trapezoidal Rule approximation
% for the integral of f(x) from x=A to x=B, using N subintervals, where
% F_STR is the string representation of f.

I=0;
g = inline(f_str);
h = (b-a)/n;

I = I + g(a);

for ii = (a+h):h:(b-h)
    I = I + 2*g(ii);
end

I = I + g(b);
I = I*h/2;

Simpson's Rule Matlab Script

function I = simprule(f_str, a, b, n)
%SIMPRULE Simpson's rule integration.
% I = SIMPRULE(F_STR, A, B, N) returns the Simpson's rule approximation
% for the integral of f(x) from x=A to x=B, using N subintervals, where
% F_STR is the string representation of f.
% An error is generated if N is not a positive, even integer.

I=0;
g = inline(f_str);
h = (b-a)/n;

if((n > 0) && (rem(n,2) == 0))
    I = I + g(a);

    for ii = (a+h):2*h:(b-h)
        I = I + 4*g(ii);
    end
   
    for kk = (a+2*h):2*h:(b-2*h)
        I = I + 2*g(kk);
    end
   
    I = I + g(b);

    I = I*h/3;
else
    disp('Incorrect Value for N')
end

Other Matlab Scripts I've Created Include:

>>> Romberg Numerical Integration - Matlab Script

>>> Lagrange Method and Newton Divided Difference Method - Matlab Scripts

>>> Newtons Method of finding Roots - Matlab Script

>>> Bisection Method of finding Roots - Matlab Script

>>> Secant Method of finding Roots - Matlab Script