CSU Vision CSU Mission Core Values CSU IGA

Transforming lives by CSU is committed to transform the lives Productivity Compassion Competent
of people and communities through high Accessibility Accountability
quality instruction and innovative Relevance
Self-disciplined
Excellence
research, development, production and

Republic of the Philippines

Cagayan State University

COLLEGE OF ENGINEERING

Carig Sur, Tuguegarao City

DEPARTMENT OF CHEMICAL ENGINEERING

Advance Engineering Mathematics for ChE

(ChE 56)

Second Semester 2016 – 2017

Course Topic: NUMERICAL INTEGRATION USING MATHLAB ®

Course Activity: RESEARCH

Name of Students: SINGSON, MARIELLA B.

TELAN, MELODY MAE B.

TULIAO, JON YOSHE S.

VENTULA, LEA B.

UBINA, KATE HYACINTH G.

Program: BSChE

Year Level: III

Date Submitted: May 12, 2017

Instructor: Engr. CAESAR P. LLAPITAN Rating:________

Date Checked: ________

1
 Table of Contents
Numerical Integration
I. Introduction ---------------------------------------------------------------------------- 1

II. Theoretical Background --------------------------------------------------------------1-3

A. Basis Engineering Integration

B. Newton Cotes Integration
a. Trapezoidal Rule
b. Simpson’s Rule
C. Gauss Quadrature
a. Romberg Integration
b. Adaptive Quadrature
c. Three Point Quadrature

III. Numerical Analysis

A. Basic Engineering Integration -------------------------------------------------------- 3-5

a. Algorithm (Pseudo-code)
b. Flowchart
c. Application
d. M-File

B. Newton Cotes Integration -------------------------------------------------------------5-17

a. Trapezoidal Rule
i. Algorithm
ii. Flowchart
iii. Application
iv. M-File
b. Simpson’s Rule
i. Algorithm
ii. Flowchart

2
 iii. Application
iv. M-File
C. Gauss Quadrature ------------------------------------------------------------------------------ 18-28
a. Romberg Integration
i. Algorithm
ii. Flowchart
iii. Application
iv. M-File
b. Adaptive Quadrature
i. Algorithm
ii. Flowchart
iii. Application
iv. M-File
c. Three-Point Quadrature
i. Algorithm
ii. Flowchart
iii. M-File

IV. Generalizations -------------------------------------------------------------------------- 29

V. References --------------------------------------------------------------------------------- 30

3
I. Introduction
Mathematically, definite integration stands for the integral of the function f (x) with
respect to the independent variable x, evaluated between the limits x =a to x =b. For functions
lying above the x axis, the integral of the function corresponds to the area under the curve of f
(x) between x = a and b. Under this is the numerical integration which is sometimes referred
as quadrature. This is an archaic term that originally meant the construction of a square
having the same area as some curvilinear figure. Today, the term quadrature is generally
taken to be synonymous with numerical definite integration.
Integration has many engineering and scientific applications. Many specific examples of
such applications could be given in all fields of engineering and science. A number of
examples relate directly to the idea of the integral as the area under a curve. For example, a
common application is to determine the mean or average of a continuous function. That mean
of a continuous function has hundreds of engineering and scientific applications. For
example, it is used to calculate the center of gravity of irregular objects in mechanical and
civil engineering and it is also used to determine the root-mean-square current in electrical
engineering. Integrals are also employed by engineers and scientists to evaluate the total
amount or quantity of a given physical variable. The integral may be evaluated over a line, an
area which is called area integral, and for the volume which is called volume integral .When
the functions to be analyzed are simple, you will normally choose to evaluate them
analytically. However, it is often difficult or impossible when the function is complicated, as
is typically the case in more realistic examples. In addition, the underlying function is often
unknown and defined only by measurement at discrete points. For both these cases,
numerical techniques are used to obtain approximate values for integrals.

II. Theoretical Background

The first numerical technique under numerical integration is basic engineering
integration. The simplest method is basic integration itself in evaluating a problem.
Applying the properties and rules of integration, a problem is solved using
b
y   f ( x)dx (1)
a

Another method discussed under this is the relationship between summation and
integration with the formula

4
 b

 f ( x)dx
Mean  a
(2)
ba
This is usually used in getting a variable given a set of values within a set of intervals.

The Newton-Cotes formulas are the most common numerical integration schemes. They are
based on the strategy of replacing a complicated function or tabulated data with a polynomial
that is easy to integrate. The trapezoidal rule is the first of the Newton-Cotes closed
integration formulas. The trapezoidal rule works by approximating the region under the graph
of the function f(x) as a trapezoid and calculating its area or simply

 f (b)  f (a)
x  a dx
b
I    f (a) 
a 
ba  (3)
Integrating the function, the equation will be simplified into
f (a)  f (b)
I  (b  a)
2 (4)
One way of improving the accuracy of the trapezoidal rule is to divide the integration interval
into a number of segments and apply the method to each segment arriving at the formula
n 1
f ( x 0 )  2 f ( x i )  f ( x n )
I  (b  a) i 1
(5)
2n
which is known to be the composite trapezoidal rule.

Aside from applying the trapezoidal rule with finer segmentation, another way to obtain a
more accurate estimate of an integral is to use higher-order polynomials to connect the points.
If there are two points equally spaced between f (a) and f (b), the four points can be
connected with a third-order polynomial. The formulas that result from taking the integrals
under these polynomials are called Simpson’s rules. The first application of Simpson’s rule
is the Simpson’s 1/3 rule that can be used when the equation to be evaluated is a second-order
with the formula.
X2 ( x  x1 )( x  x 2 ) ( x  x0 )( x  x 2 ) ( x  x0 )( x  x1 )
I  [ f ( x0 )  f ( x1 )  f ( x 2 )]dx
X0 ( x0  x1 )( x0  x 2 ) ( x1  x0 )( x1  x 2 ) ( x 2  x0 )( x 2  x1 ) (6)
simplifying this, we get
h
I  [ f ( x0 )  4 f ( x1 )  f ( x2 )] (7)
3

5
The second application of the Simpson’s rule is the Simpson’s 3/8 rule. This is similarly
derived with the 1/3 Simpson’s rule, the only difference is that this equation is used in third-
order equations with the formula
3h
I [ f ( x0 )  3 f ( x1 )  3 f ( x2 )  f ( x3 )]
8 (8)
Gauss quadrature is the name for a class of techniques to implement a strategy where
in that the constraint of fixed base points was removed and we were free to evaluate the area
under a straight line joining any two points on the curve. By positioning these points wisely,
we could define a straight line that would balance the positive and negative errors. Hence,
these cases results in a small error. The first method under gauss quadrature is the Romberg
integration with a general formula of
4 1
I I (h2 )  I (h1 )
3 3 (9)
and for computation of the error, we use

I1, k  I 2, k 1
a 
I1, k
(10)
The second application of the gauss quadrature is the two-point gauss quadrature. Just as the
trapezoidal rule, by assuming that the equation fits the integral of a constant and a linear
function exactly to fit the integral of a parabolic and a cubic function, arriving to the formula
 1   1 
I  f  f 
 3  3 (11)
Next is the three-point gauss quadrature which is a higher version of the two-point gauss
quadrature. Adaptive quadrature is another application of the gauss quadrature which shows
the desired accuracy of the answer by taking small steps in variations and larger intervals are
used when the function changes gradually.

III. Numerical Analysis

A. Basic Engineering Integration
Basic engineering Integration is simply applying the rules of integration itself in
evaluating problems like solving for mass, given a set of values and solving for the
average using summation through integration. It is the fundamental process in numerical
integration.

6
General Flowchart: Algorithm:

Start Step 1: Start

Step 2: Input 𝑥
Input 𝒙
Step 3: Use equation (1)

b Step 4: Print answer

 f ( x)dx
a Step 5: Stop

Print answer

Stop

The flowchart shows how to evaluate a definite integral given a function f(x). It has input
x, a and b.

Application:

As a specified in the following table, the earth’s density varies as a function of the distance
from its (r = 0):

r, km 0 1100 1500 2450 3400 3630

𝜌(g/cm3) 13 12.4 12 11.2 9.7 5.7

r, km 4500 5380 6060 6280 6380

𝜌(g/cm3) 5.2 4.7 3.6 3.4 3

Use numerical integration to estimate the earth’s mass (in metric tonnes) and average
density (in g/cm3). Develop vertically stacked subplots of (top) density versus radius, and
(bottom) mass versus radius. Assume that the earth is a perfect sphere.

7
Analytical Method:

6380
M  4  r (r
2
)dr
0

4r 4
M
4
M  6.1  10 24 g

M-File to implement Basis Engineering Integration:

x = [0 1100 1500 2450 3400 3630 4500 5380 6060 6280 6380];
y = [13 12.4 12 11.2 9.7 5.7 5.2 4.7 3.6 3.4 3];
rm = x*10^3;
rho = (y*10^(-3))/10^(-6);
f = 4*pi*rho.*rm.^2;
m = trapz(rm,f);
mass = max(m)

B. Newton Cotes Integration

a. Trapezoidal Rule
Trapezoidal rule is equivalent to approximating the area of the trapezoid under the
straight line connecting f (a) and f (b). To compute for the area under a curve, used the
formula width multiplied by the average height since the trapezoid is on the side.

Algorithm:

Step 1: Start

Step 2: Define and Declare function

Step 3: Input initial boundary value, final boundary value and length of interval

8
Step 4: Calculate number of strips, n = (final boundary value –final boundary
value)/length of interval.

Step 5: Perform following operation in loop

x[i]=x0+i*h

y[i]=f(x[i])

print y[i]

Step 6: Initialize se=0, s0=0

Step 7: Do the following using loop

If i %2 = 0

So=s0+y[i]

Otherwise

Se=se+y[i]

9
Flowchart:
Start

Define Function

Declare Variables

Input xo, h and xn

n  ( xn  xo ) / h

i=0

i<=n

xi   x0  i  h
yi   f ( xi )

i++

Print y[i]
se = s0 = 0
i=1

i<=n

i<=n

se= s0 + y[i] s0 = s0 + y[i]

ans  h / 3  ( y0  yn  4  s0  2  se

Print the result

Stop

10
 Trapezoidal method is based on the principle that the area under a curve which is to be
calculated is divided into a number of small segments. The bounding curve in the segment is
considered to be a straight line as a result the small enclosed area becomes a trapezium. The
area of each small trapezium is calculated and summed up. This idea is the working
mechanism of the give flowchart.

Application:

Evaluate the following integral:

4

 (1  e
x
)dx
0

i. single application of the trapezoidal rule

ii. composite trapezoidal rule with n =2 and 4

Analytical Method:

i. Trapezoidal Rule
f (a)  f (b)
b
I   [ f (a)  ( x  a)]dx
a
ba
f (a)  f (b)
I  (b  a)
2
f (0)  f (4)
I  (4  0)
2
0  0.99966
I  ( 4  0) x
2
I  1.99932

ii. Composite trapezoidal rule

n = 2, f(0) = 0, f(2) = 0.981684, f(4) = 0.99966
n 1
f ( x 0 )  2 f ( x i )  f ( x n )
I  (b  a) i 1

2n
f (0)  2 f (2)  f (4)
I  (4  0)
4
I  1.99932

n = 4 f(0) = 0, f(1) = 0.86466, f(2) = 0.98168, f(3) = 0.99752, f(4)

11
 0  2 f (1)  f (2)  f (3)  f (4)
I  (4 x)
8
I  3.34369

M-file to implement Trapezoidal rule:

function
I = trap(func,a,b,n,varargin)
% trap: composite trapezoidal rule quadrature
% I = trap(func,a,b,n,pl,p2,...):
% composite trapezoidal rule
% input:
% func = name of function to be integrated
% a, b = integration limits
% n = number of segments (default = 100)
% pl,p2,... = additional parameters used by func
% output:
% I = integral estimate
if nargin<3,error('at least 3 input arguments required'),end
if ~(b>a),error('upper bound must be greater than lower'),end
if nargin<4|isempty(n),n=100;end
x = a; h = (b - a)/n;
s=func(a,varargin{:});
for i = 1 : n-1
x = x + h;
s = s + 2*func(x,varargin{:});
end
s = s + func(b,varargin{:});
I = (b - a) * s/(2*n);

function y=f(x)
y=1-exp(-x);

12
13
Determine the distance traveled from the following velocity data:

t 1 2 3.25 4.5 6 7 8 8.5 9 10

v 5 6 5.5 7 8.5 8 6 7 7 5

Use the trapezoidal rule. In addition, determine the average velocity.

Analytical Method:

a.
56 6  5.5 5.5  7 7  8.5
d  (2  1)  (3.25  2)  (4.5  3.25)  (6.  4.5)
2 2 2 2
8.5  6 66 67 77 75
 (7  6)  (8  7)  (8.5  8)  (9.3  8.5)  (10  9.3)
2 2 2 2 2
d  58.425

M-File:

format long
t= [1 2 3.25 4.5 6 7 8 8.5 9.3 10];
v= [5 6 5.5 7 8.5 6 6 7 7 5];
p= polyfit (t,v,3)
tt = linspace(1,10);
vv = polyval(p,tt);
plot(tt,vv,t,v,'o')
xlabel('Time')
ylabel('Velocity')

14
15
b. Simpson’s Rule
Simpson’s rule is a method of numerical integration which is more accurate than the
Trapezoidal rule. It divides the area under the function to be integrated, f(x), into vertical
strips, but instead of joining the point’s f(x) with straight lines, every set of three such
successive points is fitted with a parabola.

Algorithm:

Step 1:Start (Given a function f(x))

Step 2: Get user inputs
Input
a,b=endpoints of interval
n=number of intervals (Even)

(Do the integration)

Step: Set h= (b-a)/n.

Step 4: Set sum=0.
Step 5: Begin For i= 1 to n -1

Set x =a + h*i.

If i%3=0

Then Set sum=sum+2*f(x)

Else
Set sum=sum+3*f(x) End For
Step 6: Set sum = sum + f(a)+f(b)
Step 7: Set ans = sum*(3h/8)
Set ans= sum*(h/3)
Step : End

16
 Flowchart:

Start

Enter a,b,n

Sum=0

h=(b-a)/n

Do i = 1 to n-
1

x= a +h*i

No
Is
i%2=0
?
sum=sum+3*f(x)
Yes

sum=sum+2*f(x)

I=n-1

Next

sum = sum + f(a)+f(b)

Ans = sum*(h/3)

Print ans End

17
 The flowchart shows how to implement Simpson’s 1/3 and 3/8 rule. It has inputs a,b and
n which is substituted to equations (7) and (8) to determine the area under a given curve. The
computed values will have a small percent error.

Application:

Evaluate the integral

6
1
 1  xdx
0

Using single applications of Simpson’s 1/3 rule

Analytical Method:
Equation of Simpson’s 1/3 rule

h
I  [ f ( x0 )  4 f ( x1 )  f ( x2 )]
3
n ba
n  60
n6
ba
h
n
60
h
6
h 1
X 0 1 2 3 4 5 6
Y 1 0.5 0.3333 0.2500 0.2 0.1667 0.1429

Therefore:

h
I  [ f ( x 0 )  4 f ( x1 )  f ( x 2 )]
3
h
I  [( y1  y 7 )  2( y 3  y 5 )  4( y 2  y 4  6 )]
3
I  11.75

18
M-File to implement Simpson’s Rule:

%Simpson's 3/8th and 1/8 Rule for Evaluation of Definite Integrals

using namespace std;
double f(double x)
{
double a = 1/(1+x^2); %write the function whose definite
integral is to be calcuated here
return a;
}
int main()
{ [Link](4); %set the precision
[Link](ios::fixed);
int n,i; %n is for subintervals and i is for loop
double a,b,c,h,sum=0,integral;
cout<<"\nEnter the limits of integration,\n\nInitial limit,a= ";
cin>>a;
cout<<"\nFinal limit, b="; %get the limits of
integration
cin>>b;
cout<<"\nEnter the no. of subintervals(IT IS A MULTIPLE OF 3),
\nn=; %get the no. of subintervals
cin>>n;
double x[n+1],y[n+1];
h=(b-a)/n; %get the width of the
subintervals

for (i=0;i<n+1;i++)
{ %loop to evaluate x0,...xn and y0,...yn
x[i]=a+i*h; %and store them in arrays
y[i]=f(x[i]);
}
for (i=1;i<n;i++)
{
if (i%3==0)
sum=sum+2*y[i];
else
sum=sum+3*y[i];
}
integral=3*h/8*(y[0]+y[n]+sum);
%3h/8*[y0+yn+3*(y1+y2+y4+...)+2*(y3+y6+y9+...+)]
cout<<\nThe definite integral is "<<integral<<"\n"<<endl;
return 0;
}

19
%Question: Evaluate the integral 1/(1+x) within limits 0 to 6
clc;
clear all;
close all;

f=@(x)1/(1+x);
a = 0; b = 6;
n = b-a;
h=b-a/n;
p=0;

for i=a:b
p=p+1;
x(p)=i;
y(p)=1/(1+i)
end

l=length(x);
x
y
answer=(h/3)*((y(1)+y(7))+2*(y(3)+y(5))+4*(y(2)+y(4)+y(6)))

20
C. Gauss Quadrature
a. Romberg Integration
The first application of the gauss quadrature is the Romberg integration. This technique is
designed to attain efficient numerical integrals. Through this method ad mathematical
manipulation, superior results are attained with less effort.

Flowchart:

Algorithm:
Start
Step 1: Start
Step 2: Define and declare function
Step 3: Input value to the formula
2
2 1
given to determine the I   x  
1 x
integral
Step 4: Calculate the percent error
Step 5: Determine the number of Determine the integral
iterations
Step 6: End
4 1
I I (h2 )  I (h1) A
3 3

Xi+1

I1, k  I 2 , k 1
a 
I1, k

Error<0.5% A
?

End

21
Application:
Develop an M-file function to implement Romberg integration based on Fig. 20.2. Test
the function by using it to determine the integral of 0.2+25x-200x2+675x3-900x4+400x5
from a=0 and b=0.8.

Analytical Method:
n  1, h  0.8
f (0.8)  f (0.0)
I  (0.8  0.0)  0.1728
2
 t  89.5%
n  2, h  0.4
f (0.8)  2 f (0.4)  f (0)
I  (0.8  0.0)  1.0688
4
 t  34.9%
n  4, h  0.2
f (0.8)  2 f (0.6)  2 f (0.4)  2 f (0.2)  f (0)
I  (0.8  0)  1.4848
4
 t  9.5%

Using extrapolation
4(1.0688)  (0.1728)
I  1.36747
3
 t  16.6%
4(1.4848)  (1.0688)
I  1.62347
3
 t  1.0%

Romberg Integration
4 k I (h / 2)  I (h)
I
4k 1
O(h 2 k  2 )
16(1.62347)  (1.36747)
I  1.64503
15
 t  0.0%

22
First Iteration
O( h 2 ) O( h 4 ) O( h 6 ) O(h 8 )
0.172800 1.367467
1.068800

Second Iteration
O( h 2 ) O( h 4 ) O( h 6 ) O(h 8 )
0.172800 1.367467 1.640533
1.068800 1.623467
1.484800

Third Iteration
O( h 2 ) O( h 4 ) O( h 6 ) O(h 8 )
0.172800 1.367467 1.640533 1.640533
1.068800 1.623467 1.640533
1.484800 1.639467
1.600800

M-File to implement Romberg Integration:

function I = trap(func,a,b,n,varargin)
% composite trapezoidal rule
% input:
% func = name of function to be integrated
% a, b = integration limits
% n = number of segments (default = 100)
% output:
% I = integral estimate
if nargin<3,error('at least 3 input arguments required'),end
if ~(b>a),error('upper bound must be greater than lower'),end
if nargin<4|isempty(n),n=100;end
x = a; h = (b - a)/n;
s=func(a,varargin{:});
for i = 1 : n-1
x = x + h;
s = s + 2*func(x,varargin{:});
end
s = s +func(b,varargin{:});
I = (b - a) * s/(2*n);

23
 b. Adaptive Quadrature

Adaptive Quadrature is a process in which the integral of

a function is approximated using static quadrature rules on adaptively refined
subintervals of the integration domain. Generally, adaptive algorithms are just as efficient
and effective as traditional algorithms for "well behaved" integrands, but are also
effective for "badly behaved" integrands for which traditional algorithms fail.

Flowchart:

Start
Algorithm:

Step 1: Start
Step 2: Define and declare function 2
2 1
Step 3: Input value to the formula I   x  
1 x
given to determine the
integral
Determine the integral
Step 4: Calculate the percent error
Step 5: Determine the number of
iterations 4 1
I I (h2 )  I (h1)
Step 6: End 3 3

Xi+1

I1, k  I 2 , k 1
a  End
I1, k

24
Application:

Develop an M-file function to implement adaptive quadrature based on Fig. 20.6. Test the
function by using it to determine the integral of of 0.2+25x-200x2+675x3-900x4+400x5
from a=0 and b=0.8.

Analytical Method:
n  1, h  0.8
f (0.8)  f (0.0)
I  (0.8  0.0)  0.1728
2
 t  89.5%
n  2, h  0.4
f (0.8)  2 f (0.4)  f (0)
I  (0.8  0.0)  1.0688
4
 t  34.9%
n  4, h  0.2
f (0.8)  2 f (0.6)  2 f (0.4)  2 f (0.2)  f (0)
I  (0.8  0)  1.4848
4
 t  9.5%

Using extrapolation
4(1.0688)  (0.1728)
I  1.36747
3
 t  16.6%
4(1.4848)  (1.0688)
I  1.62347
3
 t  1.0%

Romberg Integration
4 k I (h / 2)  I (h)
I
4k 1
O(h 2 k  2 )
16(1.62347)  (1.36747)
I  1.64503
15
 t  0.0%

25
M-File to implement Adaptive Quadrature:

function q = quadadapt(f,a,b,tol,varargin)
%Evaluates definite integral of f(x) from a to b
if nargin < 4 | isempty(tol),tol = 1.e-6;end
c = (a + b)/2;
fa = feval(f,a,varargin{:});
fc = feval(f,c,varargin{:});
fb = feval(f,b,varargin{:});
q = quadstep(f, a, b, tol, fa, fc, fb, varargin{:});
end
function q = quadstep(f,a,b,tol,fa,fc,fb,varargin)
% Recursive subfunction used by quadadapt.
h = b - a; c = (a + b)/2;
fd = feval(f,(a+c)/2,varargin{:});
fe = feval(f,(c+b)/2,varargin{:});
q1 = h/6 * (fa + 4*fc + fb);
q2 = h/12 * (fa + 4*fd + 2*fc + 4*fe + fb);
if abs(q2 - q1) <= tol
q = q2 + (q2 - q1)/15;
else
qa = quadstep(f, a, c,tol, fa, fd, fc, varargin{:});
qb = quadstep(f, c, b, tol, fc, fe, fb, varargin{:});
q = qa + qb;
end
end

26
c. Three Point Gauss Quadrature
Gauss-Legendre quadrature interpolates the integrand by a polynomial and integrates
the polynomial. Instead of uniformly spaced points, Gauss-Legendre uses optimally-
spaced points.

Algorithm

Step I: Start

Step II: Define integral properties

Step III: Initialize W(n,i),X(n,i)

ba 1 ba ba


b
Step IV:
a
f ( x)dx 
2  1
(
2
x
2
)dx

Step V: for i=1 to N do: Ans=ans+W(N,i)*A(X(N,i));

Step VI: Return Ans;

Step VII: End

27
Flowchart:

Start

a,b:integration Interval

N: Number of Points

f(x): FUnction Formula

Initialize W(n,i), X(n,i)

Ans=0;

ba 1 ba ba


b

2  1 2
f ( x)dx  ( x )dx
a 2

for i=1 to N do:

Ans=ans+W(N,i)*A(X(N,i));

Return Ans;

End

28
Application:

Evaluate the following integral with the three-point Gauss Quadrature formula

I   xe2 x dx
3

Analytical Method:

True Solution

3
I   xe2 x dx
0
I  504.5399

b ba 1 ba ba

a f (x) dx 
2 1
(
2
)(
2
)dx

3 1 3 3
I   ( x  ) dx
2 1 2 2
3  5   3 3  3 8   3  3 5  3 3 3   
I    f         f   0      f      
2  9   2 5  2  9   2  2  9   2 5 2   
3 5 8 5 
I   f  0.338   f (1.5)  f (2.6619) 
2 9 9 9 
3
I  (0.369  26.78  303.400)
2
I  495.8235

504.5399  495.8235
 100
495.8235
 1.73%

29
M-File to implement Three Point Quadrature:

% MATLAB script for numerically integrating f(x) over [a, b]

% using 2-point and 3-point gaussian quadrature
%
% Inputs:
% a = left endpoint of interval
% b = right endpoint of interval
% N = number of subintervals
%
% Required m-file
% f.m is m-file for function f(x)
%
% Outputs:
% G2 = value from 2-point gaussian quadrature
% G3 = value from 3-point gaussian quadrature

i=0; N=1; a=0; b=1;

% Compute exact solution

exact=exp(1)-1;

% Use MATLAB integration command

%exact=quad('f',a,b);

fprintf('\nExact answer = %e\n',exact);

30
for n=1:5
N=N*2;
i=i+1;
k=N-1;
x1=linspace(a,b,N);
h=x1(2)-x1(1);

% 2-point gaussian quadrature

sum=0;
for j=1:N-1
sum=sum+f(x1(j)+0.5*h);
end;
Gauss 2(i)=h*sum;

% 3-point gaussian quadrature

c1=0.5*h*(1-1/sqrt(3));
c2=0.5*h*(1+1/sqrt(3));
sum= 0;
for j=1:N-1
sum=sum+f(x1(j)+c1)+f(x1(j)+c2);
end;
Gauss 3(i)=0.5*h*sum;

err G2(i)=abs(G2(i)-exact); err G3(i)=abs(G3(i)-exact);

fprintf('\nNumber of subintervals = %4i \n',K(i));

fprintf('\nGauss 2 = %e Gauss 2 = %e \n',G2(i), G3(i));
fprintf('\nGauss Error 2 = %e Gauss Error 3 = %e \n\n',errG2(i),
errG3(i));

end;

loglog(K,errG2,'-or')
hold on
loglog(K,errG3,'-ob')
grid on
title('Blue = Gauss 3-point Red = Gauss 2-
point','Color','r','FontSize',14)
xlabel('Number of Subintervals','FontSize',14);
ylabel('Error','FontSize',14);

31
IV. Generalizations

Although most mathematical functions can be differentiated analytically, the same cannot
be said for integration. There are no general rules for integrating, as there are for
differentiating. We therefore need numerical methods for evaluating integrals. The common
methods used are the newton cotes formula which includes Trapezoidal rule and Simpson’s
1/3 and 3/8 rule, and Gauss Quadrature which includes Romberg Integration, Adaptive
Quadrature, Two Point Quadrature, and Three Point Quadrature. These methods are used to
estimate the area under a curve. For trapezoidal rule, its error is large compared to Simpson’s
rule. Therefore Simpson’s 1/3 rule is more accurate than trapezoidal rule. Simpson’s 1/3 rule
is third order accurate even though it is based on only three points. In other words, it yields
exact results for cubic polynomials even though it is derived from a parabola. In trapezoidal
rule, if the function being integrated is linear, it will be exact because the second derivative of
a straight line is zero. Otherwise, for functions with second- and higher-order derivatives
(with curvature), some error can occur.

When the function is known and high accuracy is required, a method such as Romberg
integration occurs. For Romberg integration, that technique can be used to generate an
integral estimate within a pre specified error tolerance. It is a method for combining two
numerical integral estimates to obtain a third, more accurate value. Gauss-quadrature
formulas employ x values that are positioned between the integration limits in such a manner
that a much more accurate integral estimate results. Although Romberg integration is more
efficient than the composite Simpson’s 1/3 rule, both use equally spaced points. This
constraint does not take into account that some functions have regions of relatively abrupt
changes where more refined spacing might be required. Hence, to achieve a desired accuracy,
fine spacing must be applied everywhere even though it is only needed for the regions of
sharp change. Adaptive quadrature methods remedy this situation by automatically adjusting
the step size so that small steps are taken in regions of sharp variation.

32
V. References

Chapman, S. J. (2001). MATLAB Programming for Engineers.

Chapra, S. C. (2005). Applied Numerical Methods with MATLAB for Engineers and
Scientist. New York: McGraw-Hill Icn.

Hahn, B. D., & Valentine, D. T. (2001). Essential MATLAB for Engineers and Scientist.
Oxford: Elsvier Ltd.

33