/* * Program 4.4 Gauss_QD.C * Purpose: Demonstrate how to find the Integration of given function * by Using Gauss-Legendre Quadrature. * Author : Wachara Rodsumrid. * Office : Dept of Physics. Rajamangala Institute of Technology. * Language: Turbo C v.2.0 * Last update : 6 November 95 **************************************************************/ #include #include double YourFunc(double); double GetLowerLimit(void); double GetUpperLimit(void); int GetNumberofNode(void); int GaussQuadrature(int,double,double,double*); void IntegrationResult(int,double*); void main() { double *result; double lower_limit, upper_limit; int n=0 /* number of point to be used for integrate */; int error; clrscr(); printf("\t\tFinding the Integration using Gauss Quadrature.\n"); printf ("\t\t==============================================\n"); lower_limit = GetLowerLimit(); upper_limit = GetUpperLimit(); n = GetNumberofNode(); error = GaussQuadrature(n,lower_limit,upper_limit,result); IntegrationResult(error,result); printf("\n\nYour work was done. Goodbye.\n"); } double GetLowerLimit() { char buffer[255]; printf("\n\tInput the Lower Limit of Integration : "); gets(buffer); return(atof(buffer)); } double GetUpperLimit() { char buffer[255]; printf("\n\tInput the Upper Limit of Integration : "); gets(buffer); return(atof(buffer)); } int GetNumberofNode() { char buffer[255]; printf("\n\tInput the number of Gauss point to be used (from 2 to 6) : "); gets(buffer); return(atoi(buffer)); } void IntegrationResult(int err,double *result) { int i; switch (err) { case 0 : { printf("\nThe Result of Integration using Gauss Quadrature = %f",*result); break; } case 1 : { printf("\nThe number of node is not between 2 to 6.\n"); printf("or Lower limit equals to upper limit.\n"); break; } } /* switch (err) */ } /*================== Function Start here. ==========================*/ int GaussQuadrature(int n,double lower, double upper,double *integration) /* Function Description: Compute the integral of given function from lower limit to upper limit using the Gauss Quadrature method. User's Define Variable : None User's Define Function : double YourFunc( double ); Input : n -- number of node for integration from 1 to 6 lower -- lower limit. upper -- upper limit. integration -- variable for keeping the result Output : the integration of f(x) in integration. Return : 0 -- if no error. 1 -- if 6 < n < 2 or lower limit = upper limit *************************************************************************/ { double c[7] ; /* Coefficient of function */ double x[7],u[7] ; /* node point of function */ double fvalue; int i,j; if (n < 2 || n > 6 ) return 1; if ((upper - lower) == 0) return 1; switch (n) { case 1 : break; case 2 : { u[2] = 0.5773502691; c[1] = c[2] = 1; u[1] = -u[2]; break; } case 3 : { u[2] = 0; u[3] = 0.7745966692; u[1] = -u[3]; c[2] = 0.8888888888; c[1] = c[3] = 0.5555555555; break; } case 4 : { u[3] = 0.3399810435; u[2] = -u[3]; u[4] = 0.8611363115; u[1] = -u[4]; c[3] = 0.6521451548; c[4] = 0.3478548451; c[1] = c[4]; c[2] = c[3]; break; } case 5 : { u[3] = 0; u[4] = 0.5384693101; u[5] = 0.9061798459; u[1] = -u[5]; u[2] = -u[4]; c[3] = 0.5688888888; c[2] = c[4] = 0.4786286704; c[1] = c[5] = 0.2369268850; break; } case 6 : { u[4] = 0.2386191860; u[5] = 0.6612093864; u[6] = 0.9324695142; u[1] = -u[6]; u[2] = - u[5]; u[3] = -u[4]; c[4] = 0.4679139345; c[5] = 0.3607615730; c[6] = 0.1713244923; c[1] = c[6]; c[2] = c[5]; c[3] = c[4]; break; } } /* switch (n) */ /* ******* Used for studying or debuging program ******* printf( "------------------------------------------------------------------ \n" ); printf( " Term Coef node x f(x) \n" ); printf( "------------------------------------------------------------------ \n" ); ********************************************************/ for( i = 1; i <= n; i++ ){ x[i] = 0.5*( (upper + lower) + (upper - lower)*u[i]); } *integration = 0; for( i = 1; i <= n; i++ ){ fvalue = YourFunc(x[i]); /* Calculate functional values */ *integration += fvalue*c[i]; /* printf( " %3d %12.5e %12.5e %12.5e %12.5e\n",i,c[i],u[i], x[i], fvalue); */ } *integration *= (upper - lower)/2; /* Don't forget dx = (b-a)du/2 */ return 0; } /* =============== The End of Gauss Quadrature =============*/ double YourFunc(double x) /* Write your function which you want to intergrate in here */ { return(2*sqrt(1-(x-1)*(x-1))); } /*=======================================================*/