Docsity
Docsity

Prepare for your exams
Prepare for your exams

Study with the several resources on Docsity


Earn points to download
Earn points to download

Earn points by helping other students or get them with a premium plan


Guidelines and tips
Guidelines and tips

Numerical Integration, Exams of Integrated Case Studies

Numerical Integration with so easy and usefull doc fotr university students

Typology: Exams

2016/2017

Uploaded on 11/19/2017

oguz-biyiklioglu
oguz-biyiklioglu 🇹🇷

1 document

1 / 6

Toggle sidebar

This page cannot be seen from the preview

Don't miss anything!

bg1
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)+2*sum(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-2
% Prentice-Hall Pub. Inc.
% One Lake Street
% Upper Saddle River, NJ 07458
h=(b-a)/(2*M);
1
pf3
pf4
pf5

Partial preview of the text

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;