# Problems on Numerical Integration – C Programming

In my previous posts, I showed you guys how to write C programs for various Numerical Integration Techniques, like Trapezoidal Rule, and Simpson’s 1/3 & 3/8 Rule.

I have also written quite a few posts on C Programs for Numerical Root Finding techniques.

So, in this post we will be solving some problems based on the above knowledge, and thus it will be a good exercise to write some complex programs applying various Numerical Techniques together.

### Exercise 1:

Integrate to an accuracy of 1 in 10^5 for given limits a and b:
$\int_5^{10} \frac{\arctan x}{x^2}dx$
Use Trapezoidal & Simpson rule.

### PROGRAM:

/**********************************
********PROBLEM 6.4.1**************
**********************************/
#include<stdio.h>
#include<math.h>
double f(double x){
return atan(x)/(x*x);
}
/*Function definition to perform integration by Trapezoidal Rule */
double trapezoidal(double f(double x), double a, double b, int n){
double x,h,sum=0,integral;
int i;
h=fabs(b-a)/n;
for(i=1;i<n;i++){
x=a+i*h;
sum=sum+f(x);
}
integral=(h/2)*(f(a)+f(b)+2*sum);
return integral;
}
/*Function definition to perform integration by Simpson's 1/3rd Rule */
double simpsons(double f(double x), double a,double b,double n){
double h,integral,x,sum=0;
int i;
h=fabs(b-a)/n;
for(i=1;i<n;i++){
x=a+i*h;
if(i%2==0){
sum=sum+2*f(x);
}
else{
sum=sum+4*f(x);
}
}
integral=(h/3)*(f(a)+f(b)+sum);
return integral;
}
main(){
int n,i=2;
double a,b,h,x,integral,eps,integral_new;

/*Ask the user for necessary input */
printf("\nEnter the initial limit: ");
scanf("%lf",&a);
printf("\nEnter the final limit: ");
scanf("%lf",&b);
printf("\nEnter the desired accuracy: ");
scanf("%lf",&eps);
integral_new=simpsons(f,a,b,i);

/* Perform integration by simpson's 1/3rd for different number of sub-intervals until they converge to the given accuracy:*/
do{
integral=integral_new;
i=i+2;
integral_new=simpsons(f,a,b,i);
}while(fabs(integral_new-integral)>=eps);

printf("\nThe integral using Simpson's Rule is: %lf for %d sub-intervals.\n",integral_new,i);

i=2;
/* Perform integration by trapezoidal rule for different number of sub-intervals until they converge to the given accuracy:*/
do{
integral=integral_new;
i++;
integral_new=trapezoidal(f,a,b,i);
}while(fabs(integral_new-integral)>=eps);

printf("The integral using Trapezoidal Rule is: %lf\n with %d intervals",integral_new,i);
}



### EXERCISE 2:

The time period of a pendulum is given by the integral
$T=4\int_0^{\pi/2} \frac{1}{1-\sin^2(A/2) \sin^2x}dx$
where A is the amplitude of oscillations. For small amplitudes it is possible to approximate the
time period to
$T_1=2\pi\left[ 1 + \left( \frac{A}{4} \right)^2\right]$
Plot T, T1 and the percentage difference between T and T1 as functions of A for $0\leq A\leq \pi$.

### PROGRAM:

/***************************************
***********PROBLEM 6.4.2****************
***************************************/
#include<stdio.h>
#include<math.h>
/*Function definition to perform integration by Simpson's 1/3rd Rule */
double simpsons(double A,double f(double x,double A),double a,double b,double n){
double h,integral,x,sum=0;
int i;
h=fabs(b-a)/n;
for(i=1;i<n;i++){
x=a+i*h;
if(i%2==0){
sum=sum+2*f(x,A);
}
else{
sum=sum+4*f(x,A);
}
}
integral=(h/3)*(f(a,A)+f(b,A)+sum);
return integral;
}
double f(double x, double A){
return 1/(1-sin(A/2)*sin(A/2)*sin(x)*sin(x));
}
double T1(double A){
return 2*M_PI*(1+pow(A/4,2));
}
double T(double f(double x,double A),double A,int n){
double integral;
integral=simpsons(A,f,0,M_PI/2,n);
return 4*integral;
}
main(){
FILE *fp=NULL;
fp=fopen("integrationProblem2.txt","w");
double A;
double t,t1,diffT;
for(A=0;A<=M_PI;A=A+0.1){
t=T(f,A,600);
t1=T1(A);
diffT=(t-t1)/t*100;
fprintf(fp,"%lf\t%lf\t%lf\t%lf\n",A,t,t1,diffT);
}
}


### GnuPlot Command:

plot './integrationProblem2.txt' w l t "T", '' u 1:3 w l t "T1", '' u 1:4 w l t "% Error"

### EXERCISE 3:

Locate the smallest positive root of the function F(x) , given by:
$F(x)=\int_0^\pi\cos \left[x^a\cos(t)\right]\sin^{2n+1}tdt$
to an accuracy of 4 decimal places, for n = 1 and a = 1.5.

### PROGRAM:

/***************************************
***********PROBLEM 6.4.4****************
***************************************/
#include<stdio.h>
#include<math.h>
/*Function definition to perform integration by Simpson's 1/3rd Rule */
double simpsons(double x,double f(double t,double x),double a,double b,int n){
double h,integral,t,sum=0;
int i;
h=fabs(b-a)/n;
for(i=1;i<n;i++){
t=a+i*h;
if(i%2==0){
sum=sum+2*f(t,x);
}
else{
sum=sum+4*f(t,x);
}
}
integral=(h/3)*(f(t,a)+f(t,b)+sum);
return integral;
}
double f(double t,double x){
return cos(pow(x,1.5)*cos(t))*pow(sin(t),3);
}
double F(double x){
return simpsons(x,f,0,M_PI,600);
}
/*The following function performs the bisection procedure and also prints the values of various variables at each iteration*/
double printBisection(double f(double x),double a, double b, double eps, int  maxSteps){
double c;
if(f(a)*f(b)<=0){
int iter=1;
/*Bisection Method begins that tabulates the various values at each iteration*/
printf("____________________________________________________________________________________\n");
printf("iter\ta\t\tb\t\tc\t\tf(c)\t\t|a-b|\n");
printf("____________________________________________________________________________________\n");
do{
c=(a+b)/2;
printf("%d.\t%lf\t%lf\t%lf\t%lf\t%lf\n",iter,a,b,c,f(c),fabs(a-b));
if(f(a)*f(c)>0){
a=c;
}
else if(f(a)*f(c)<0){
b=c;
}
iter++;

}while(fabs(a-b)>=eps&&iter<=maxSteps);
printf("___________________________________________________________________________________________________\n");
return c;
}
else{
printf("\nSorry! Either the root doesn't exist in the given interval or there are multiple roots in this interval.\nPlease enter a different set of guesses.\n");
return 9999;
}
}
main(){
//Let us first tabulate the function for a given range of x
double xmin, xmax;
printf("Enter the lower value for x:\nxmin = ");
scanf("%lf",&xmin);
printf("Enter the upper value for x:\nxmax = ");
scanf("%lf",&xmax);
double x;
printf("x\t\tf(x)\n");
printf("__________________________\n");
for(x=xmin;x<=xmax;x=x+0.1){
printf("%lf\t%lf\n",x,F(x));
}
char choice='y';
while(choice=='y'){
//Begin Bisection Routine
printf("Begining Bisection Routine:\n");
double a,b,eps;
int maxSteps;
printf("Enter the initial guess:\na = ");
scanf("%lf",&a);
printf("b = ");
scanf("%lf",&b);
printf("Enter the desired accuracy:");
scanf("%lf",&eps);
printf("Enter the maximum no. of iterations to be performed: ");
scanf("%d",&maxSteps);
double root=printBisection(F,a,b,eps,maxSteps);
//9999 is the error code returned by the bisection function if the given interval dosen't bracket the root or contains more than 1 root
if(root!=9999){
printf("One of the roots of the function in the given interval is: %lf",root);
}

printf("\nDo you want to find more roots?\ny/n\n");
scanf(" %c", &choice);
}
}



### EXERCISE 4:

Use the integral representation of the Bessel function:
$J_n(z)=\frac{1}{2\pi}\int_0^{2\pi}\cos(z\cos x))dx$
to find its zeroes in the range $0\leq z \leq12$ by secant method.

### PROGRAM:

/***************************************
***********PROBLEM 6.4.5****************
***************************************/
#include<stdio.h>
#include<math.h>
/*Function definition to perform integration by Simpson's 1/3rd Rule */
double simpsons(double z,double f(double x,double z),double a,double b,int n){
double h,integral,x,sum=0;
int i;
h=fabs(b-a)/n;
for(i=1;i<n;i++){
x=a+i*h;
if(i%2==0){
sum=sum+2*f(x,z);
}
else{
sum=sum+4*f(x,z);
}
}
integral=(h/3)*(f(a,z)+f(b,z)+sum);
return integral;
}
double f(double x,double z){
return cos(z*cos(x));
}
double J0(double z){
return 1/(2*M_PI)*simpsons(z,f,0,2*M_PI,1000);
}
/*Secant Method Function that tabulates the values at each iteration*/
double printSecant(double f(double x), double x1, double x2, double eps, int maxSteps){
int iter=1;
double x3;
printf("___________________________________________________________________\n");
printf("iter\tx1\t\tx2\t\tx3\t\tf(x3)\n");
printf("___________________________________________________________________\n");
do{
x3=(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1));
printf("%d\t%lf\t%lf\t%lf\t%lf\n",iter,x1,x2,x3,f(x3));
x1=x2;
x2=x3;
iter++;
}while(fabs(f(x3))>eps&&iter<=maxSteps);
printf("___________________________________________________________________\n");
return x3;
}
main(){
double zmin,zmax,z;
//Let's tabulate the function in a given range
printf("Enter the range of the z values between which the zeroes are to be determined:\nzmin = ");
scanf("%lf",&zmin);
printf("zmax = ");
scanf("%lf",&zmax);
printf("z\t\tJ0(z)\n");
printf("_______________________________\n");
for(z=zmin;z<=zmax;z=z+0.1){
printf("%lf\t%lf\n",z,J0(z));
}
//Begin Secant-Routine of Root Finding
double a,b,eps;
int maxSteps;
printf("**Secant Method Root finding routine begins**\n");
char choice='y';
while(choice=='y'){
printf("Enter the initial guesses for the Secant Method:\na = ");
scanf("%lf",&a);
printf("b = ");
scanf("%lf",&b);
printf("Enter the accuray desired:");
scanf("%lf",&eps);
printf("Enter the maximum no. of iterations:");
scanf("%d",&maxSteps);
double root=printSecant(J0,a,b,eps,maxSteps);
printf("One of the roots is %lf\n",root);
printf("Do you want to find more roots?y/n\n");
scanf(" %c", &choice);
}
}


### EXERCISE 5:

The spherical Bessel function of order n is given by
$j_n(z)=\frac{z^n}{2^{n+1}n!}\int_0^\pi \cos (z\cos \theta)\sin^{2n+1}\theta d\theta$
Find all the roots of j2(z) between 0 and 10.

### PROGRAM:

/***************************************
***********PROBLEM 6.4.6****************
***************************************/
#include<stdio.h>
#include<math.h>
//Function to perform factorial
double factorial(int n){
int i;
double fact=1;
for(i=n;i>=1;i--){
fact=fact*i;
}
return fact;
}
/*Function definition to perform integration by Simpson's 1/3rd Rule */
double simpsons(double z,int n,double f(double theta,double z,int n),double a,double b,int intervals){
double h,integral,theta,sum=0;
int i;
h=fabs(b-a)/intervals;
for(i=1;i<intervals;i++){
theta=a+i*h;
if(i%2==0){
sum=sum+2*f(theta,z,n);
}
else{
sum=sum+4*f(theta,z,n);
}
}
integral=(h/3)*(f(a,z,n)+f(b,z,n)+sum);
return integral;
}
double f(double theta,double z, int n){
return cos(z*cos(theta))*pow(sin(theta),2*n+1);
}
double Jn(double z, int n){
return pow(z,n)/(pow(2,n+1)*factorial(n))*simpsons(z,n,f,0,M_PI,1000);
}
double J2(double z){
return Jn(z,2);
}
/*Secant Method Function that tabulates the values at each iteration*/
double printSecant(double f(double x), double x1, double x2, double eps, int maxSteps){
int iter=1;
double x3;
printf("___________________________________________________________________\n");
printf("iter\tx1\t\tx2\t\tx3\t\tf(x3)\n");
printf("___________________________________________________________________\n");
do{
x3=(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1));
printf("%d\t%lf\t%lf\t%lf\t%lf\n",iter,x1,x2,x3,f(x3));
x1=x2;
x2=x3;
iter++;
}while(fabs(f(x3))>eps&&iter<=maxSteps);
printf("___________________________________________________________________\n");
return x3;
}
main(){
double z,zmin,zmax,z1,h,eps,zinc=0.1;
int maxSteps;
printf("Enter the interval in which you want to find the roots of J2(z):\nzmin = ");
scanf("%lf",&zmin);
printf("zmax = ");
scanf("%lf",&zmax);
printf("Enter the step-width for tabulation: ");
scanf("%lf",&h);
printf("The values of J2(z) in the given range are:\n");
printf("z\t\tJ2(z)\n");
printf("______________________________________\n");
for(z=zmin;z<=zmax;z=z+h){
printf("%lf\t%lf\n",z,J2(z));
}
printf("Enter the accuracy desired for the roots: ");
scanf("%lf",&eps);
printf("Enter the max. no. of iterations that are to be performed: ");
scanf("%d",&maxSteps);
printf("The roots in the given range are: \n");
for(z=zmin;z<=zmax;z=z+zinc){
z1=z+zinc;
if(J2(z)*J2(z1)<0){
printf("\nIn the interval: %lf and %lf\n",z,z1);
double root=printSecant(J2,z,z1,eps,maxSteps);
printf("Root: %lf\n",root);
}
}
}


### REFERENCES:

The above problems have been taken from the Computer Programming & Numerical Analysis Manual by Dr. Shobhit Mahajan.

[wpedon id="7041" align="center"]

## 5 thoughts on “Problems on Numerical Integration – C Programming”

1. can you help to write c++ program for interpolation and extrapolation

2. would you please provide me with a soft copy of Computer Programming & Numerical Analysis Manual

3. can u explain about M_PI, because the varibale is undefine

This site uses Akismet to reduce spam. Learn how your comment data is processed.