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.

PhD researcher at Friedrich-Schiller University Jena, Germany. I'm a physicist specializing in theoretical, computational and experimental condensed matter physics. I like to develop Physics related apps and softwares from time to time. Can code in most of the popular languages. Like to share my knowledge in Physics and applications using this Blog and a YouTube channel.



3 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 *