/******************************************************** * Integration by Gauss method of a real function F=F(X) * * or F=F(X,Y) or F=F(X,Y,Z). The integral is calculated * * by using from 2 to 10 Gauss points. * * ----------------------------------------------------- * * Ref.: "Mécanique des vibrations linéaires By * * M. Lalanne, P. Berthier and J. Der Hagopian, * * Masson, Paris, 1980" [BIBLI 16]. * * ----------------------------------------------------- * * SAMPLE RUN: * * * * (Integrate F=sin(x) from x=0 to x=1). * * * * INTEGRATION OF A REAL FUNCTION BY GAUSS * * F(X), F(X,Y) or F(X,Y,Z) * * * * Number of variables (1 to 3): 1 * * * * How many Gauss points (2 to 10): 4 * * * * Minimum value of X: 0 * * Maximum value of X: 1 * * * * Value of integral = 0.45969769 * * * * C++ version by J-P Moreau, Paris. * * (www.jpmoreau.fr) * ********************************************************/ #include #include // Labels: e10,e20,e30,e40 double A[11], H[11]; double x,x1,x2,y,yy1,y2,z,z1,z2; int i,j,k,n,n1,n2; double a1,a2,a3,b1,b2,b3,xi9; //define here function F double F(double x,double y,double z) { return sin(x); } void main() { printf("\n INTEGRATION OF A REAL FUNCTION BY GAUSS\n"); printf(" F(X), F(X,Y) or F(X,Y,Z)\n"); printf("\n Number of variables (1 to 3): "); scanf("%d",&n2); printf("\n How many Gauss points (2 to 10): "); scanf("%d",&n); printf("\n"); n1 = n - 1; switch(n1) { case 1:A[1] = -0.57735026919; H[1] = 1.0; break; case 2:A[1] = -0.774596669241; A[2] = 0.0; H[1] = 0.555555555556; H[2] = 0.888888888889; break; case 3:A[1] = -0.8611363115939999; A[2] = -0.339981043585; H[1] = 0.347854845137; H[2] = 0.652145154863; break; case 4:A[1] = -0.906179845939; A[2] = -0.538469310106; A[3] = 1e-12; H[1] = 0.236926885056; H[2] = 0.478628670499; H[3] = 0.568888888889; break; case 5:A[1] = -0.9324695142029999; A[2] = -0.661209386466; A[3] = -0.238619186083; H[1] = 0.171324492379; H[2] = 0.360761573048; H[3] = 0.467913934573; break; case 6:A[1] = -0.949107912343; A[2] = -0.741531185599; A[3] = -0.405845151377; A[4] = 0.0; H[1] = 0.129484966169; H[2] = 0.279705391489; H[3] = 0.381830050505; H[4] = 0.417959183673; break; case 7:A[1] = -0.949107912343; A[2] = -0.741531185599; A[3] = -0.405845151377; A[4] = 0.0; H[1] = 0.129484966169; H[2] = 0.279705391489; H[3] = 0.381830050505; H[4] = 0.417959183673; break; case 8:A[1] = -0.960289856497; A[2] = -0.796666477414; A[3] = -0.525532409916; A[4] = -0.183434642496; H[1] = 0.10122853629; H[2] = 0.222381034453; H[3] = 0.313706645878; H[4] = 0.362683783378; break; case 9:A[1] = -0.968160239508; A[2] = -0.836031107327; A[3] = -0.6133714327000001; A[4] = -0.324253423404; A[5] = 0.0; H[1] = 0.0812743883616; H[2] = 0.180648160695; H[3] = 0.260610696403; H[4] = 0.31234707704; H[5] = 0.330239355001; break; case 10:A[1] = -0.973906528517; A[2] = -0.865063366689; A[3] = -0.679409568299; A[4] = -0.433395394129; A[5] = -0.148874338982; H[1] = 0.0666713443087; H[2] = 0.149451349151; H[3] = 0.219086362516; H[4] = 0.26926671931; H[5] = 0.295524224715; } for (i = 1; i<=n/2; i++) { j = n + 1 - i; A[j] = -A[i]; H[j] = H[i]; } printf("\n Minimum value of X: "); scanf("%lf", &x1); printf("\n Maximum value of X: "); scanf("%lf", &x2); printf("\n"); if (n2 - 1 > 0) { printf("\n Minimum value of Y: "); scanf("%lf", &yy1); printf("\n Maximum value of Y: "); scanf("%lf", &y2); printf("\n"); } if (n2 - 2 > 0) { printf("\n Minimum value of Z: "); scanf("%lf", &z1); printf("\n Maximum value of Z: "); scanf("%lf", &z2); } xi9 = 0; for (i = 1; i<=n; i++) { if (n2 - 1 == 0) goto e10; for (j = 1; j<=n; j++) { if (n2 - 2 == 0) goto e10; for (k=1; k<=n; k++) { e10: a1 = -(x1 - x2) / 2; b1 = (x1 + x2) / 2; x = a1 * A[i] + b1; if (n2 - 1 == 0) goto e20; a2 = -(yy1 - y2) / 2; b2 = (yy1 + y2) / 2; y = a2 * A[j] + b2; if (n2 - 2 == 0) goto e20; a3 = -(z1 - z2) / 2; b3 = (z1 + z2) / 2; z = a3 * A[j] + b3; e20: switch(n2) { case 1:xi9 += F(x,y,z) * H[i] * a1; goto e40; break; case 2:xi9 += F(x,y,z) * H[i] * H[j] * a1 * a2; goto e30; break; case 3:xi9 += F(x,y,z) * H[i] * H[j] * H[k] * a1 * a2 * a3; } } e30: ;} e40:;} printf("\n Value of integral = %12.8f\n\n", xi9); } // end of file tgauss.cpp