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 funi-e(x) x."3-x.*cos (x) ; fun2- ( .*exp (-2.*x.*2) ; Number of steps in each methods N 100; Limit of integration for function al-0; bl-2*pi; %Limit of integration for function 2. a2--l; b2-1; fprintf(n nFor 1st function f (x)-' disp (fun1) %Exact integration for function 1 I_ext - integral (funl, al,bl); fprint f ( "\nExact integration value for a=%f to b-%f 1s %f. n',al,bl,I ext) %Left Riemann integration for function 1 I_rmn-left_riemann (funl,al,bl,N); fprintf( "\n\tLeft Riemann integration for function 1 for N=%d is %f. n',N,Irmn) tprint f ( "Error in Left Riemann integration is %e.\n',abs(1-ext-Irmn )) VTrapí zoidal integration for function 1 I_trp-trapizoidal( funl,al,bl,N) tprint f( "\n\ tTrapizoidal integration for function 1 for N=2d is %f. n,N,I_trp) fprint f( "Error in Trapizoidal integration is %e,\n', abs (1-ext-1-trp)) Simpson integration for function 1 I smp-simpson( funl,al, b1,N)i tprint f ( .\n\ tSimpson integration for function 1 for N-%d is %f. n,N, I_smp) fprint f ( Error in Simpson integration is %e.\n', abs (1-ext-1-smp)) Romberg integration for function 1 I_rom-romberg ( fun1,al,bl,10); fprint f ( \n\tRomberg integration for function 1 for N-ed is %f. n',10,I rom) fprint f ( "Error in Romberg integration is %e.\n', abs (1-ext-1-rom)) fprintf(n nFor 2nd function f (x)-' disp (fun2) %Exact integration for function 2 I-ext = integral (fun2, a2,b2); fprint f ( "\nExact integration value for a=%f to b-%f 1s %f. n',a2,b2,I ext) %Left Riemann integration for function 2

I rmn-left riemann(fun2, a2,b2,N); fprint f ( .\n\ tLeft Riemann integration for function 2 for N-%d is %f. n',N, I rmn) tprint f (Error in Left Riemann integration is %e.\n.abs(1-ext-1-rmn)) Trapizoidal integration for function 2 Itrp-trapizoidal ( fun2,a2, b2,N) fprint f ( .\n\tTrapizoidal integration for function 2 for N-%d is %f. n,N,Itrp) fprint f ( "Error in Trapizoidal integration is %e.\n", abs (1-ext-1-trp)) Simpson integration for function 2 I smp-simpson(fun2, a2, b2, N)i fprint f ( .\n\tSimpson integration for function 2 for N-%d 1s n',N,Ismp) fprint f ( Error in Simpson integration is %e.\n', abs (1-ext-1-smp)) Romberg integration for function 2 I_rom-romberg(fun2,a2, b2,10); fprint f ( "\n\tRomberg Integration for function 2 for N=2d is %f n',10,I rom) fprint f ("Error in Romberg integration is %e.\n', abs (1-ext-1-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-l:length (xx)-1 val-val+dx double (func (xx1)); end end %%Matlab function for Trapizoldal 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

xx 1 =xx (1); val-val+dx*double (func (xxl)); end val-val+dx* (0.5*double (func (xx(1)))+0.5*double (func (xx (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 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) fMatlab inline function ⅝ a, b- integration interval % n - number of rows in Romberg tableau to 11- to 11 erence of 2 succesive romberg value % Output: r - Romberg tableau containing the computed values of the integral % h - step length % n - size of Romberg table r(1,1) - (b- a) * (func(a) + func(b)) / 2; for j -2:N subtotal -0; for i- 1:2* (j-2) subtotalsubtotal +func (a (2 * i - 1) * h(j)); end r(j, i) = 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)i end val-r(end, end);

용옹옹%%옹옹%%옹옹%%옹옹%%옹옹 End of Code 용융 응용 For 1st function f (x)- e (x)x. 3-x. *cos (x) Exact integration value for a-0.000000 to b-6.283185 İS 389.63636. Left Riemann integration for function 1 for N-100 is 382.079993. Error in Left Riemann integration is 7.556372e+00. Trapizoidal integration for function1 for N-100 is 389.675328. Error in Trapizoidal integration is 3.896364e-02. Simpson integration for function 1 for N-100 is 389.636364 Error in Simpson integration is 0.000000e+00. Romberg integration for function 1 for N-10 is 389.636364. rror in Romberg integration is 5.684342e-14 For 2nd function f(x)- (x)4. *exp(-2.*x. *2) Exact integration value for a=-1.000000 to b=1.000 OOO İS 4.785 152. Left Riemann integration for function 2 for N-100 is 4.785008 Error in Left Riemann integration is 1.443538e-04 Trapizoidal integration for function 2 for N-100 is 4.785008. Error in Trapizoidal integration is 1.443538-04 Simpson integration for function 2 for N-100 is 4.785152 Error in Simpson integration is 1.542456e-08. Romberg integration for function 2 for N-10 is 4.785152. Error in Romberg integration is 0.000000e+00. Published with MATLAB R2018a

%%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 %%%%%%%%%%%%%%%%