Numerical Integration Formulas

Here are the M-files to implement composite trapezoidal rule for equally spaced data and unequally spaced data. 

Composite Trapezoidal Rule (Equally Spaced Data) M-file

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);

Composite Trapezoidal Rule (Unequally Spaced Data) M-file

function I = trapuneq(x,y)

% trapuneq: unequal spaced trapezoidal rule quadrature

% I = trapuneq(x,y):

% Applies the trapezoidal rule to determine the integral

% for n data points (x, y) where x and y must be of the

% same length and x must be monotonically ascending

% input:

% x = vector of independent variables

% y = vector of dependent variables

% output:

% I = integral estimate

if nargin<2,error('at least 2 input arguments required'),end

if any(diff(x)<0),error('x not monotonically ascending'),end

n = length(x);

if length(y)~=n,error('x and y must be same length'); end

s = 0;

for k = 1:n-1

s = s + (x(k+l)-x(k))*(y(k)+y(k+l))/2;

end

I = s;

Numerical Integration Formulas

Problem 19.13.

The total mass of a variable density rod is given by 

where m = mass, ρ(x) = density, A_c(x) = cross-sectional area, x = distance along the rod, and L = total length of rod. The following data have been measured for a 20-m (L) length rod. Determine the mass in grams to the best possible accuracy. 

Solution:

First, be aware of the units. Convert the unit for distance into cm to ensure the calculation is dimensionally homogenous. Also, compute the product of density and cross-sectional area of the rod. 

To attain best possible accuracy for the integral, we can apply combination of trapezoidal rule (1st segment), Simpson's 1/3 rule (2nd & 3rd segments) and Simpson's 3/8 rule (last three segments). 

Numerical Integration Formulas

Problem 19.5.

The function f(x) = e^-x can be used to generate the following table of unequally spaced data: 

Evaluate the integral from a = 0 to b = 1.2 using (a) analytical means, (b) the trapezoidal rule, and (c) a combination of the trapezoidal and Simpson's rules wherever possible to attain the highest accuracy. For (b) and (c), compute the true percent relative error. 

Solution:

(a) The analytical solution can be evaluated as

(b) Apply trapezoidal rule to each segment and sum the results

(c) Apply trapezoidal rule to the first segment (h = 0.1), Simpson's 3/8 rule to the next three segments (h = 0.2), and Simpson's 1/3 rule to the last two segments (h = 0.25). Then, sum the results

Polynomial Interpolation: Inverse Interpolation

Problem 17.9.

Employ inverse interpolation to determine the value of x that corresponds to f(x) = 0.93 for the following tabulated data: 

Note that the values in the table were generated with the function f(x) = x²/(1+x²).

(a) Determine the correct value analytically. 

(b) Use quadratic interpolation and the quadratic formula to determine the value numerically. 

(c) Use cubic interpolation and bisection to determine the value numerically. 

Solution:

(a) Using the function available, manipulate it to determine corresponding x value:

(b) A quadratic interpolating polynomial can be fit to the last three points: (3, 0.9), (4, 0.941176), and (5, 0.961538) using polyfit function:

>> format long

>> x = [3 4 5];

>> fx = [0.9 0.941176 0.961538];

>> a = polyfit(x,fx,2)

a =

  -0.010407000000000   0.114024999999999   0.651588000000002

Thus, the best-fit parabola that passes through the last three points is f₂(x) = -0.010407x² + 0.114025x + 0.651588. 

We must therefore find the root of 

0.93 = -0.010407x² + 0.114025x + 0.651588

or

0 = -0.010407x² + 0.114025x - 0.278412

The quadratic formula can be used to calculate:

Thus, the second root, 3.67295, is a good estimate of the true value. 

(c) A cubic interpolating polynomial can be fit to the last four points: (2, 0.8), (3, 0.9), (4, 0.941176), and (5, 0.961538) using polyfit function:

>> format long

>> x = [2 3 4 5];

>> fx = [0.8 0.9 0.941176 0.961538];

>> a = polyfit(x,fx,3)

a =

   0.006335000000000  -0.086427000000002   0.411770000000007   0.271487999999994

Thus, the best-fit cubic that passes through the last four points is f₃(x) = 0.006335x³ - 0.086427x² + 0.41177x + 0.271488. 

We must therefore find the root of 

0.93 = 0.006335x³ - 0.086427x² + 0.41177x + 0.271488

or

0 = 0.006335x³ - 0.086427x² + 0.41177x - 0.658512

Bisection method (bisect.m) can be employed to find the root of this polynomial. Using initial guesses of x_l = 3 and x_u = 4 (that bracket the root), a value of 3.61884 is obtained.

>> fx = @(x) 0.006335*x^3 - 0.086427*x^2 + 0.41177*x - 0.658512;

>> [root,fx,ea,iter]=bisect(fx,3,4)

root =

   3.618841171264648

fx =

     1.082876788238707e-08

ea =

     5.270606093347705e-05

iter =

    19

Polynomial Interpolation

Problem 17.5.

Given the data

Calculate f(4) using Newton's interpolating polynomials of order 1 through 4. Choose your base points to attain good accuracy. That is, the points should be centered around and as close as possible to the unknown. What do your results indicate regarding the order of the polynomial used to generate the data in the table? 

Solution:

The points can be ordered so that they are close to and centered around the unknown. A divided-difference table can then be developed as:

Note that the fact that the fourth divided difference is zero means that the data was generated with a third-order polynomial.

Polynomial Interpolation

Problem 17.4.

Given the data 

(a) Calculate f(3.4) using Newton's interpolating polynomials of order 1 through 3. Choose the sequence of the points for your estimates to attain the best possible accuracy. That is, the points should be centered around and as close as possible to the unknown. 

(b) Repeat (a) but use the Lagrange polynomial. 

Solution:

(a) Newton's interpolating polynomial. Ordering of points:

Note that based purely on the distance from the unknown (nearer to 3.4), the fourth point would be (2,5). However, because it provides better balance and is located only a little bit farther from the unknown, the point at (5,0) is chosen.

Set up divided difference table to ease the evaluation of each order of polynomial interpolation. 

First divided differences:

Second divided differences:

Third divided difference:

Thus, the divided difference table is

First order (linear) interpolation:

Second order (quadratic) interpolation:

Third order (cubic) interpolation:

(b) Lagrange interpolating polynomial.

First order interpolation:

Second order interpolation:

Third order interpolation:

Polynomial Interpolation

Here are the M-files to implement Newton interpolation and Lagrange interpolation.

Newton Interpolation M-file

function yint = Newtint(x,y,xx)
% Newtint: Newton interpolating polynomial
% yint = Newtint(x,y,xx): Uses an (n - 1)-order Newton
% interpolating polynomial based on n data points (x, y)
% to determine a value of the dependent variable (yint)
% at a given value of the independent variable, xx.
% input:
% x = independent variable
% y = dependent variable
% xx = value of independent variable at which
% interpolation is calculated
% output:
% yint = interpolated value of dependent variable

% compute the finite divided differences in the form of a
% difference table
n = length(x);
if length(y)~=n, error('x and y must be same length'); end
b = zeros(n,n);
% assign dependent variables to the first column of b.
b(:,1) = y(:); % the (:) ensures that y is a column vector.
for j = 2:n
for i = 1:n-j+1
b(i,j) = (b(i+1,j-1)-b(i,j-1))/(x(i+j-1)-x(i));
end
end
% use the finite divided differences to interpolate
xt = 1;
yint = b(1,1);
for j = 1:n-1
xt = xt*(xx-x(j));
yint = yint+b(1,j+1)*xt;
end

Lagrange Interpolation M-file

% Lagrange interpolating polynomial based on n data points

% to determine a value of the dependent variable (yint) at

% a given value of the independent variable, xx.

% input:

% x = independent variable

% y = dependent variable

% xx = value of independent variable at which the

% interpolation is calculated

% output:

% yint = interpolated value of dependent variable

n = length(x);

if length(y)~=n, error('x and y must be same length'); end

s = 0;

for i = 1:n

product = y(i);

for j = 1:n

if i ~= j

product = product*(xx-x(j))/(x(i)-x(j));

end

end

s = s+product;

end

yint = s;