% ### EXintegration1.m ### 2016.01.13 CB % Numerical integration example % (modified from) Original source: % http://ef.engr.utk.edu/ef230-2011-01/modules/matlab-integration/ clear; % ---------------------- % User parameters F = @(x)(sin(x)); % function to integrate %F = @(x)(exp(-x.^2/2)); % function to integrate xL= [0 pi]; % integration limits N= 2; % Method A - # of "rectangles" for LEFT and RIGHT pts= [3 4 5 6 10 25]; % Method B - # of points to consider integrating (via trapz function) dur= 1; % Method B - pause duration [s] for trapz loop % ---------------------- % *************** % Show the curve figure(1); fplot(F,[xL(1),xL(2)]) % a quick way to plot a function xlabel('x'); ylabel('F(x)'); % *************** % Method A % Approximate the integral via brute force LEFT and RIGHT Riemann sums sumL= 0; sumR=0; delX= (xL(2)-xL(1))/N; % step-size x= linspace(xL(1),xL(2),N+1); % add one since N is # of 'boxes' and is really N-1 for nn=1:N sumL= sumL + F(x(nn))*delX; sumR= sumR + F(x(nn+1))*delX; end disp(['left-hand rule yields =',num2str(sumL),' (for ',num2str(N+1),' steps)']); disp(sprintf('right-hand rule yields = %g', sumR)); disp(['area calculated by 1/2*(LEFT+RIGHT) for ',num2str(N),' points =',num2str(0.5*(sumL+sumR))]); % *************** % Method B % Approximate the integral via trapz for different numbers of points for np=pts figure(2); clf % clear the current figure hold on % allow stuff to be added to this plot x = linspace(xL(1),xL(2),np); % generate x values y = F(x); % generate y values a2 = trapz(x,y); % use trapz to integrate % Generate and display the trapezoids used by trapz for ii=1:length(x)-1 px=[x(ii) x(ii+1) x(ii+1) x(ii)]; py=[0 0 y(ii+1) y(ii)]; fill(px,py,ii) end fplot(F,[xL(1),xL(2)]); xlabel('x'); ylabel('F(x)'); disp(['area calculated by trapz.m for ',num2str(np),' points =',num2str(a2)]); title(['area calculated by trapz.m for ',num2str(np),' points =',num2str(a2)]); pause(dur); % wait a bit end % *************** % Method C a1 = quad(F,xL(1),xL(2)); % use quad to integrate msg = ['area calculated by quad.m = ' num2str(a1,10)]; disp(msg);