Download Numerical Integration and more Exams Integrated Case Studies in PDF only on Docsity!
Numerical Integration
We have three main ways to do numerical integration using equally spaced points: The Trapezoidal Rule, The Simpson’s Rule, and Romberg Integration. Here are examples of MatLab m-files for the Trape- zoidal Rules (traprulefh.m), Simpson’s Rule, Romberg Integration(rombfh.m), and one that illustrates the Romberg table of “improved” values (rombdemofh.m):
function I=traprulefh(Fun,a,b,n); %% function I=traprulefh(Fun,a,b,n); %% %% COMPOSITE TRAPEZOIDAL RULE FOR INTEGRATION %% %% INPUT: %% %% Fun - the function must be a function handle %% OR as a row vector of data [f(x0),f(x1),...,f(xn)] %% a - lower endpoint of integration %% b - upper endpoint of integration %% n - an integer BIGGER THAN 2 (x_n is the endpoint b, altogether %% n+1 points counting endpoints) %% OUTPUT: %% I - the approximate value of the integral %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% h=(b-a)/n; points=a+h[0:n]; Fvalues=Fun(points); I=Fvalues(1)+Fvalues(n+1)+2sum(Fvalues(2:n)); I=h*I/2; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s=simprl(f,a,b,M)
%Input - f is the integrand % - a and b are upper and lower limits of integration % - M is the number of subintervals %Output - s is the simpson rule sum
%If f is defined as an M-file function use the @ notation % call s=simprl(@f,a,b,M). %If f is defined as an anonymous function use the % call s=simprl(f,a,b,M).
% NUMERICAL METHODS: Matlab Programs % (c) 2004 by John H. Mathews and Kurtis D. Fink % Complementary Software to accompany the textbook: % NUMERICAL METHODS: Using Matlab, Fourth Edition % ISBN: 0-13-065248- % Prentice-Hall Pub. Inc. % One Lake Street % Upper Saddle River, NJ 07458
h=(b-a)/(2*M);
s1=0; s2=0;
for k=1:M x=a+h(2k-1); s1=s1+f(x); end for k=1:(M-1) x=a+h2k; s2=s2+f(x); end
s=h(f(a)+f(b)+4s1+2*s2)/3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function app=rombfh(f,a,b,n); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% function app=rombfh(f,a,b,n); %% ROMBERG INTEGRATION %% INPUT: %% f - the function which must be a function handle %% a - lower endpoint of integration %% b - upper endpoint of integration %% n - ‘‘depth’’ number of rows in the romberg table %% %% OUTPUT: %% app - approximate integral %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Step 1 h=b-a; R(1,1)=h(f(a)+f(b))/2; %Step 3 for i=2:n; p=2^(i-2); x=a-.5h+h(1:p); R(2,1)=(R(1,1)+hsum(f(x)))/2; for j=2:i; R(2,j)=R(2,j-1)+(R(2,j-1)-R(1,j-1))/(4^(j-1)-1); end; h=h/2; R(1,1:i)=R(2,1:i); end; app=R(1,n);
function R=rombdemofh(Fun,a,b,n); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% function R=rombdemofh(Fun,a,b,n); %% ROMBERG INTEGRATION %% DEMO VERSION %% INPUT: %% Fun - the function which must be entered as a function handle
simpson_result=left_simpson+right_simpson; end;
Here is an m-file for the function f (x) = (3x^2 − 1)^2 /^3
function z=thirdrootfn(x); %% function defining (3x.^2-1).^(2/3) z=(real((3x.^2-1).^(1/3))).^2;
This function does not have a derivative at x = 1/
3 in [0, 1]. We apply the adaptive Simpson’s Rule twice here and ilustrate it using the m-file
function [simpson_result,achoice]=recur_simpson_demo(Fun,a,b,TOL,level,level_max,achoice); % function [simpson_result,achoice]=recur_simpson_demo(Fun,a,b,TOL,level,level_max,achoice); % Collecting the left hand end points % Adaptive Simpson’s Rulea for the integral of Fun on [a,b] % following Cheney-Kincaid Chapter 6. % INPUT % Fun - a function m-file name in ’ ’ % a,b - lower and upper limits of integration respectively % TOL - desired accuracy of the integral % level - level of recursion (start with level=0 ); % level_max - maximum depth of recursion % OUTPUT % simpson_result - value of the integral % achoice - left end points with duplication %echo on achoice=[achoice,a]; level=level+1; h=b-a; c=(b+a)/2; one_simpson=h(feval(Fun,a)+4feval(Fun,c)+feval(Fun,b))/6; d=(a+c)/2; e=(c+b)/2; two_simpson=h(feval(Fun,a)+4feval(Fun,d)+2feval(Fun,c)+4feval(Fun,e)+feval(Fun,b))/12; if level >= level_max simpson_result=two_simpson; message=’Level Limit Exceeded’ elseif abs(two_simpson-one_simpson)<15*TOL; simpson_result=two_simpson+(two_simpson-one_simpson)/15; else [left_simpson,achoice]=recur_simpson_demo(Fun,a,c,TOL/2,level,level_max,achoice); [right_simpson,achoice]=recur_simpson_demo(Fun,c,b,TOL/2,level,level_max,achoice); simpson_result=left_simpson+right_simpson; end;
The above m-file keeps track of the endpoints used and so we see how the rule adapts. Here are three runs.
[simpson_result,achoice]=recur_simpson_demo(’thirdrootfn’,0,1,10^(-5),0,5,[]); simpson_result = 0. [simpson_result,achoice]=recur_simpson_demo(’thirdrootfn’,0,1,10^(-8),0,8,[]);
(^00) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
1
(^00) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
1
Figure 1: Recursive Simpson
x=0:.001:1; plot(x,f(x),achoice,f(achoice),’o’) simpson_result = 0. [simpson_result,achoice]=recur_simpson_demo(’thirdrootfn’,0,1,10^(-15),0,12,[]); simpson_result = 0.
The m-files given above for the adaptive Simpson’s Rule require the function to be given separately as an m-file. Here are programs modified for use with function handles:
function simpson_result=recur_simpsonfh(Fun,a,b,TOL,level,level_max); %%function simpson_result=recur_simpsonfh(Fun,a,b,TOL,level,level_max); % Adaptive Simpson’s Rulea for the integral of Fun on [a,b] % following Cheney-Kincaid Chapter 6. % INPUT % Fun - a function as a function handle % a,b - lower and upper limits of integration respectively % TOL - desired accuracy of the integral % level - level of recursion % level_max - maximum depth of recursion % OUTPUT % simpson_result - value of the integral %echo on level=level+1; h=b-a; c=(b+a)/2; one_simpson=h(Fun(a)+4Fun(c)+ Fun(b))/6; d=(a+c)/2; e=(c+b)/2; two_simpson=h(Fun(a)+4Fun(d)+2Fun(c)+4Fun(e)+Fun(b))/12; if level >= level_max simpson_result=two_simpson; message=’Level Limit Exceeded’ elseif abs(two_simpson-one_simpson)<15*TOL; simpson_result=two_simpson+(two_simpson-one_simpson)/15; else left_simpson=recur_simpsonfh(Fun,a,c,TOL/2,level,level_max); right_simpson=recur_simpsonfh(Fun,c,b,TOL/2,level,level_max); simpson_result=left_simpson+right_simpson; end;
