Question

These instructions are written with the assumption that code will be done in matlab. You might find the following built in commands useful: length, plot, xlabel, ylabel, title, legend, fzero, plot, disp, axis, axes, min, max.

2. Numerical Integration (Quadrature). Write FOUR of your own numerical integration routines. One should use left-end, right-end, or mid-points, another should use the trapezoid method, another should use Simpson’s method, and the fourth should use either Guassian Quadrature or Romberg’s method. Use your four routines to approximate two integrals numerically so that they are accurate to 1 part in 1 billion (if possible, if this is not possible explain why). First use your routines to approximate f(x) = x 3 − x cos(x) on the interval [0, 2π], you can and should check that your answer is correct. Then use your routines to approximate f(x) = 4e −2x 2 on the interval [−1, 1]. Use one of matlab’s built in integration functions to approximate both integrals also. Each integral should be: numerically approximated FOUR ways, checked vs one of matlab’s built-in numerical integration methods, and the polynomial/trig function should be checked exactly.

Answer 1

%%Matlab code for finding integration using 4 different method
clear all
close all
%functions for which integration have to do
fun1=@(x) x.^3-x.*cos(x) ;
fun2=@(x) 4.*exp(-2.*x.^2) ;
%Number of steps in each methods
N=100;
%Limit of integration for function 1.
a1=0; b1=2*pi;
%Limit of integration for function 2.
a2=-1; b2=1;

fprintf('\n\nFor 1st function f(x)= ')
disp(fun1)
%Exact integration for function 1
I_ext = integral(fun1,a1,b1);
fprintf('\nExact integration value for a=%f to b=%f is %f.\n',a1,b1,I_ext)

%Left Riemann integration for function 1
I_rmn=left_riemann(fun1,a1,b1,N);
fprintf('\n\tLeft Riemann integration for function 1 for N=%d is %f.\n',N,I_rmn)
fprintf('Error in Left Riemann integration is %e.\n',abs(I_ext-I_rmn))

%Trapizoidal integration for function 1
I_trp=trapizoidal(fun1,a1,b1,N);
fprintf('\n\tTrapizoidal integration for function 1 for N=%d is %f.\n',N,I_trp)
fprintf('Error in Trapizoidal integration is %e.\n',abs(I_ext-I_trp))

%Simpson integration for function 1
I_smp=simpson(fun1,a1,b1,N);
fprintf('\n\tSimpson integration for function 1 for N=%d is %f.\n',N,I_smp)
fprintf('Error in Simpson integration is %e.\n',abs(I_ext-I_smp))

%Romberg integration for function 1
I_rom=romberg(fun1,a1,b1,10);
fprintf('\n\tRomberg integration for function 1 for N=%d is %f.\n',10,I_rom)
fprintf('Error in Romberg integration is %e.\n',abs(I_ext-I_rom))

fprintf('\n\nFor 2nd function f(x)= ')
disp(fun2)
%Exact integration for function 2
I_ext = integral(fun2,a2,b2);
fprintf('\nExact integration value for a=%f to b=%f is %f.\n',a2,b2,I_ext)

%Left Riemann integration for function 2
I_rmn=left_riemann(fun2,a2,b2,N);
fprintf('\n\tLeft Riemann integration for function 2 for N=%d is %f.\n',N,I_rmn)
fprintf('Error in Left Riemann integration is %e.\n',abs(I_ext-I_rmn))

%Trapizoidal integration for function 2
I_trp=trapizoidal(fun2,a2,b2,N);
fprintf('\n\tTrapizoidal integration for function 2 for N=%d is %f.\n',N,I_trp)
fprintf('Error in Trapizoidal integration is %e.\n',abs(I_ext-I_trp))

%Simpson integration for function 2
I_smp=simpson(fun2,a2,b2,N);
fprintf('\n\tSimpson integration for function 2 for N=%d is %f.\n',N,I_smp)
fprintf('Error in Simpson integration is %e.\n',abs(I_ext-I_smp))

%Romberg integration for function 2
I_rom=romberg(fun2,a2,b2,10);
fprintf('\n\tRomberg integration for function 2 for N=%d is %f.\n',10,I_rom)
fprintf('Error in Romberg integration is %e.\n',abs(I_ext-I_rom))

%%Matlab function for Left Riemann integration
function val=left_riemann(func,a,b,N)
    % func is the function for integration
    % a is the lower limit of integration
    % b is the upper limit of integration
    % N number of rectangles to be used
    val=0;
    %splits interval a to b into N+1 subintervals
    xx=linspace(a,b,N+1);
    dx=xx(2)-xx(1); %x interval
    %loop for Riemann integration
        for i=1:length(xx)-1
            xx1=xx(i);
            val=val+dx*double(func(xx1));
        end
end

%%Matlab function for Trapizoidal integration
function val=trapizoidal(func,a,b,N)
    % func is the function for integration
    % a is the lower limit of integration
    % b is the upper limit of integration
    % N number of rectangles to be used
    val=0;
    %splits interval a to b into N+1 subintervals
    xx=linspace(a,b,N+1);
    dx=xx(2)-xx(1); %x interval
    %loop for Riemann integration
        for i=2:length(xx)-1
            xx1=xx(i);
            val=val+dx*double(func(xx1));
        end
       val=val+dx*(0.5*double(func(xx(1)))+0.5*double(func(xx(end))));
end

%%Matlab function for Simpson integration
function val=simpson(func,a,b,N)
    % func is the function for integration
    % a is the lower limit of integration
    % b is the upper limit of integration
    % N number of rectangles to be used
    %splits interval a to b into N+1 subintervals
    xx=linspace(a,b,N+1);
    dx=xx(2)-xx(1); %x interval
    val=(dx/3)*(double(func(xx(1)))+double(func(xx(end))));
    %loop for Riemann integration
        for i=2:length(xx)-1
            xx1=xx(i);
            if mod(i,2)==0
                val=val+(dx/3)*4*double(func(xx1));
            else
                val=val+(dx/3)*2*double(func(xx1));
              
            end
        end    
end

%%Matlab function for Romberg integration
function val = romberg(func,a,b,N)
% f - Matlab inline function
% a,b - integration interval
% n - number of rows in Romberg tableau
% toll- tollerence of 2 succesive romberg value
% Output:
% r - Romberg tableau containing the computed values of the integral
% h - step length
% n - size of Romberg table
        h = (b - a) ./ (2.^(0:N-1));
        r(1,1) = (b - a) * (func(a) + func(b)) / 2;
        for j = 2:N
            subtotal = 0;
            for i = 1:2^(j-2)
                subtotal = subtotal + func(a + (2 * i - 1) * h(j));
            end
            r(j,1) = r(j-1,1) / 2 + h(j) * subtotal;
            for k = 2:j
                r(j,k) = (4^(k-1) * r(j,k-1) - r(j-1,k-1)) / (4^(k-1) - 1);
            end
        end
        val=r(end,end);
end

%%%%%%%%%%%%%%%%%%% End of Code %%%%%%%%%%%%%%%%