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.

SOLUTION:

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);
  
 	/*Print the answer */
  	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);

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

OUTPUT:

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 .

SOLUTION:

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"

OUTPUT:

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.

SOLUTION:

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


OUTPUT:

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.

SOLUTION:

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

OUTPUT:

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.

SOLUTION:

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

OUTPUT:

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. would you please provide me with a soft copy of Computer Programming & Numerical Analysis Manual

Leave a Reply

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

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