Numerical Differentiation and Integration with MATLAB I

Composite Trapezoidal Rule. To approximate the integral

$\displaystyle \int_a^b\approx\frac{h}{2}(f(a)+f(b))+h\sum_{k=1}^{M-1}f(x_k)
$

by sampling $ f(x)$ at the $ M+1$ equally spaced points $ x_k=a+kh$, for $ k=0,1,2,\ldots,M$. Notice that $ x_0=a$ and $ x_M=b$.
function s=traprl(f,a,b,M)
%Input    - f is the integrand input as a string 'f'
%         - a and b are upper and lower limits of integration
%         - M is the number of subintervals
%Output   - s is the trapezoidal rule sum

h=(b-a)/M;
s=0;
for k=1:(M-1)
   x=a+h*k;
   s=s+feval(f,x);
end
s=h*(feval(f,a)+feval(f,b))/2+h*s;
You are given the function $ f(x)=2+sin(2\sqrt x)$ for the interval $ [1 ,6 ]$.
  1. Plot the function.
  2. Use the composite trapezoidal rule with 11 sample points to compute an approximation to the integral of $ f(x)$ taken over $ [1 ,6 ]$ by using the MATLAB program given above.
  3. Do the error analysis. Error term for the composite trapezoidal rule is given as;

    $\displaystyle E(f,h)=-\frac{(b-a)}{12}h^2f''(\xi)=O(h^2)
$

  4. Calculate the exact value of the integration by using MATLAB. Compare your results for the aspects of integration and error analysis.
  5. Repeat the procedure with increased number of sample points.
    Solution:
    function week9lsg(fs,a,b,vec)
    % input fs as string 
    f=inline(fs); % inline
    format short;
    N=length(vec);
    disp('Exact      Approximated       Approximated-Exact      Maxerror        Minerr');
    for i=1:N
    approximated=num2str(traprl(f,a,b,vec(i)));
    fx2 = diff(fs,2);
    h=(b-a)/vec(i);
    fs = '2  + sin (2*(x^.5))';
    errmax=num2str(((b-a)/12)*h^2*fx2(a));
    errmin=num2str(((b-a)/12)*h^2*fx2(b));
    exact=int(fs,a,b);
    D=[exact,approximated,approximated-exact,errmax,errmin];
    disp(D);
    end
    
    save with the name week10lsg.m. Then;
    >> vec=[10 20 40 50 100 500 1000]';
    >> week10lsg('2  + sin (2*(x^.5))',1,6,vec);
    Exact  Approximated  Approximated-Exact   Maxerror        Minerr
    8.1835, 8.1939, .10421e-1, 1.87500, 1.5625
    8.1835, 8.1860, .25208e-2, 0.46875, 0.39062
    8.1835, 8.1841, .62079e-3, 0.11719, 0.097656
    8.1835, 8.1839, .42079e-3, 0.07500, 0.0625
    8.1835, 8.1836, .12079e-3, 0.01875, 0.015625
    8.1835, 8.1835, .20792e-4, 0.00075, 0.000625
    8.1835, 8.1835, .20792e-4, 0.0001875, 0.00015625
    


2004-12-28