Trapezoidal Rule – Algorithm, Implementation in C With Solved Examples


Numerical Methods & Algorithms / Sunday, October 21st, 2018

At first we deduce the general integration formula based on Newton’s forward interpolation formula and after that we will use it to formulate Trapezoidal Rule and Simpson’s 1/3 rd rule.

The Newton’s forward interpolation formula for the equi-spaced points xi, i =0, 1, …, n, xi = x0 + ih is

\[\phi \left( x \right)={{y}_{0}}+u\Delta {{y}_{0}}+\frac{u\left( u-1 \right)}{2!}{{\Delta }^{2}}{{y}_{0}}+\frac{u\left( u-1 \right)\left( u-2 \right)}{3!}{{\Delta }^{3}}{{y}_{0}}+….\]

\[where~~u=\frac{x-{{x}_{0}}}{h},~~h~~is~~the~~spacing.\]

Let the interval [a, b] be divided into n equal subintervals such that a = x0 < x1 < x2 < … < xn = b. Then

\[I=\int_{a}^{b}{f(x)dx=\int_{{{x}_{0}}}^{{{x}_{n}}}{\phi \left( x \right)}}dx\]

\[=\int_{{{x}_{0}}}^{{{x}_{n}}}{\left[ {{y}_{0}}+u\Delta {{y}_{0}}+\frac{u\left( u-1 \right)}{2!}{{\Delta }^{2}}{{y}_{0}}+\frac{u\left( u-1 \right)\left( u-2 \right)}{3!}{{\Delta }^{3}}{{y}_{0}}+…. \right]}dx\]

Since x = x0 + uh, dx = hdu, when x = x0 then u = 0 and x = xn then u = n. Thus,

\[I=\int_{0}^{n}{\left[ {{y}_{0}}+u\Delta {{y}_{0}}+\frac{{{u}^{2}}-u}{2!}{{\Delta }^{2}}{{y}_{0}}+\frac{{{u}^{3}}-3{{u}^{2}}+2u}{3!}{{\Delta }^{3}}{{y}_{0}}+…. \right]hdu}\]

\[=h\left[ {{y}_{0}}\left[ u \right]_{0}^{n}+\Delta {{y}_{0}}\left[ \frac{{{u}^{2}}}{2} \right]_{0}^{n}+\frac{{{\Delta }^{2}}{{y}_{0}}}{2!}\left[ \frac{{{u}^{3}}}{3}-\frac{{{u}^{2}}}{2} \right]_{0}^{n}+\frac{{{\Delta }^{3}}{{y}_{0}}}{3!}\left[ \frac{{{u}^{4}}}{4}-{{u}^{3}}+{{u}^{2}} \right]_{0}^{n}+…. \right]\]

\[=nh\left[ {{y}_{0}}+\frac{n}{2}\Delta {{y}_{0}}+\frac{2{{n}^{2}}-3n}{12}{{\Delta }^{2}}{{y}_{0}}+\frac{{{n}^{3}}-4{{n}^{2}}+4n}{24}{{\Delta }^{3}}{{y}_{0}}+…. \right]……..(1)\]

From this formula, one can generate different integration formulae by substituting n = 1, 2, 3, …

Trapezoidal Rule

Substituting n = 1 in the equation (1). In this case all differences higher than the first difference become zero. Then

\[\int_{{{x}_{0}}}^{{{x}_{n}}}{f\left( x \right)dx=h\left[ {{y}_{0}}+\frac{1}{2}\Delta {{y}_{0}} \right]}\]

\[\Rightarrow \int_{{{x}_{0}}}^{{{x}_{n}}}{f\left( x \right)dx=h}\left[ {{y}_{0}}+\frac{1}{2}\left( {{y}_{1}}-{{y}_{0}} \right) \right]=\frac{h}{2}\left( {{y}_{0}}+{{y}_{1}} \right)……….(2)\]

The formula (2) is known as the Trapezoidal Rule.

In this formula, the interval [a, b] is considered as a single interval, and it gives a very rough answer. But, if the interval [a, b] is divided into several subintervals and this formula is applied to each of these subintervals then a better approximate result may be obtained.

This formula is known as composite formula, deduced below.

Composite Trapezoidal Rule

Let the interval [a, b] be divided into n equal subintervals by the points a = x0 < x1 < x2 < … < xn = b, where xi = x0 + ih, i = 1, 2, 3, …,n.

Applying the trapezoidal rule to each of the subintervals, one can find the composite formula as

\[\int_{a}^{b}{f\left( x \right)dx=\int_{{{x}_{0}}}^{{{x}_{1}}}{f\left( x \right)dx}}+\int_{{{x}_{1}}}^{{{x}_{2}}}{f\left( x \right)dx}+….+\int_{{{x}_{n-1}}}^{{{x}_{n}}}{f\left( x \right)dx}\]

\[=\frac{h}{2}\left[ {{y}_{0}}+{{y}_{1}} \right]+\frac{h}{2}\left[ {{y}_{1}}+{{y}_{2}} \right]+….+\frac{h}{2}\left[ {{y}_{n-1}}+{{y}_{n}} \right]\]

\[\therefore \int_{a}^{b}{f\left( x \right)dx=}\frac{h}{2}\left[ {{y}_{0}}+2\left( {{y}_{1}}+{{y}_{2}}+….+{{y}_{n-1}} \right)+{{y}_{n}} \right]\]

Algorithm of Trapezoidal Rule

Step 1. Input f(x);
Step 2. Rread a,b,n; //the lower and upper limits and number of subintervals
Step 3. Compute h=(b-a)/n;
Step 4. Set sum =[f(a)+f(a+nh)]/2;
Step 5. for i=1 to n-1 do
			Compute sum = sum + f(a+ih);
		endfor;
Step 6. Compute result = sum * h;
Step 7. Print result;

Trapezoidal Rule Implementation in C

/* This program finds the value of integration of a function
   by Trapezoidal rule. Here we assume that f(x) = x^3. */

#include<stdio.h>

void main()
{
	float a,b,h,sum;
	int n,i;
	float f(float);
	printf("Enter the values of a, b: ");
	scanf("%f%f",&a,&b);
	printf("Enter the value of n: ");
	scanf("%d",&n);
	h=(b-a)/n;
	sum=(f(a)+f(a+n*h))/2;
	for(i=1;i<n;i++)
		sum+=f(a+i*h);
	sum=sum*h;
	printf("The value of the integration is %8.5f: ",sum);
}

float f(float x)
{
	return(x*x*x);
}

Output

Enter the values of a, b: 0 1
Enter the value of n: 100
The value of the integration is 0.25002

 Example 01

Find the value of

\[\int_{1.2}^{1.6}{\left( x+\frac{1}{x} \right)dx}\]

Taking 4 subintervals, correct up to four significant figures.

Solution:

\[Let~~f\left( x \right)=\left( x+\frac{1}{x} \right)\]

\[Here~~{{x}_{0}}=1.2,~~{{x}_{n}}=1.6,~~n=4\]

\[\therefore h=\frac{1.6-1.2}{4}=0.1\]

The tabulated values of f(x) for different values of x are given below:

x 1.2 1.3 1.4 1.5 1.6
f(x) 2.033333 2.069231 2.114286 2.166667 2.225

By Trapezoidal Rule, we have

\[\int_{a}^{b}{f\left( x \right)dx=}\frac{h}{2}\left[ {{y}_{0}}+2\left( {{y}_{1}}+{{y}_{2}}+{{y}_{3}} \right)+{{y}_{4}} \right]\]

\[\int_{1.2}^{1.6}{\left( x+\frac{1}{x} \right)dx}=\frac{0.1}{2}\left[ 2.033333+2\left( 2.069231+2.114286+2.166667 \right)+2.225 \right]\]

\[\therefore \int_{1.2}^{1.6}{\left( x+\frac{1}{x} \right)dx}=0.8477\]

Leave a Reply

Your email address will not be published. Required fields are marked *