/* ----------------------- MODULE kubnec.cpp ------------------------- * * "Numerical Algorithms with C, By Gisela Engeln-Muellges * * and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * ----------------------------------------------------------------------*/ #include // for REAL, ZERO, FABS, ONE, TWO, POW #include static REAL K4KnotGew (int, int); static void Kub4RoST (REAL a, REAL b, int Nx, REAL c, REAL d, int Ny, int nST, REAL f (REAL,REAL), REAL* W); static int RoRiExtr (REAL* RoFo, int nRoFo, int Ordnung, REAL* Wert, REAL* FehlerSch); static void Kub4BuST (REAL a, REAL b, int Nx, REAL c, REAL d, int Ny, int j, REAL f (REAL,REAL), REAL* L); static int BuRiExtr (REAL* BuFo, int nBuFo, int Ordnung, REAL* Wert, REAL* FehlerSch); static int BuNenner (int j); static int Kub3Nec3n (REAL Px, REAL Py, REAL Qx, REAL Qy, REAL Rx, REAL Ry, int n, REAL f(REAL,REAL), REAL* W); int Kub4NeCn ( REAL a, REAL b, int Nx, REAL c, REAL d, int Ny, int Verfahren, REAL f (REAL,REAL), REAL* Wert, int Schaetzen, REAL* FehlerSch ); /*********************************************************************** * Cubature over rectangles using Newton-Cotes formulas. * * * * Integrate the function f(x,y) over a rectangle (a,b) x (c,d) using * * the summed Newton-Cotes cubature formulas for sub-rectangles. * * * * Parameters: * * REAL a,b left, right x-end points * * int Nx number of x-intervals * * REAL c,d ditto for y * * int Ny * * int Verfahren Number of method : * * 1: trapezoidal rule * * 2: Simpson's rule * * 3: 3/8 * * 4: 4/90 * * 5: 5/288 * * 6: 6/840 * * 7: 7/17280 formula * * double f (x,y) function * * REAL *Wert value for integral * * int Schaetzen 0: no estimation * * 1: estimate error by repeating cubature for half * * the step size * * REAL *FehlerSch error estimate for Wert * * * * Return value : * * 0: o.k. * * 1: Nx improper * * 2: Ny improper * * 3: method number incorrect * * 4: Integration interval of length zero * * * * Author: Uli Eggermann, 3.31.1996 * ***********************************************************************/ /********************************************************************** * Global constants and variables * **********************************************************************/ static REAL const K_1 [] = { 1./2, 1./2 }; static REAL const K_2 [] = { 1./6, 4./6, 1./6 }; static REAL const K_3 [] = { 1./8, 3./8, 3./8, 1./8 }; static REAL const K_4 [] = { 7./90, 32./90, 12./90, 32./90, 7./90 }; static REAL const K_5 [] = { 19./288, 75./288, 50./288, 50./288, 75./288, 19./288 }; static REAL const K_6 [] = { 41./840, 216./840, 27./840, 272./840, 27./840, 216./840, 41./840 }; static REAL const K_7 [] = { 751./17280, 3577./17280, 1323./17280, 2989./17280, 2989./17280, 1323./17280, 3577./17280, 751./17280 }; static REAL const *KubArr [] = { K_1, K_2, K_3, K_4, K_5, K_6, K_7 }; int KubVer = -1; /* global variables: method number */ int KubX, KubY; /* number of x-, y-intervals */ /********************************************************************** * Cubature over rectangles via Newton-Cotes * **********************************************************************/ int Kub4NeCn (REAL a, REAL b, int Nx, REAL c, REAL d, int Ny, int Verfahren, REAL f (REAL,REAL), REAL* Wert, int Schaetzen, REAL* FehlerSch) { int i, j, k, Ordnung; REAL Hab, Hcd, Wert1 = ZERO; if (Nx < 1) return (1); if (Ny < 1) return (2); if (Verfahren < 1 || Verfahren > 7) return (3); if (a == b || c == d) return (4); /* method: 1,2,3,4,5,6,7 */ Ordnung = (Verfahren / 2) * 2 + 2; /* order: 2,4,4,6,6,8,8 */ KubVer = Verfahren; /* copy to global variable */ for (k = 1; k <= (Schaetzen ? 2 : 1); k++) { KubX = k * Nx * KubVer; KubY = k * Ny * KubVer; *Wert = ZERO; /* initialize */ Hab = (b - a) / KubX; /* step sizes */ Hcd = (d - c) / KubY; for (i = 0; i <= KubX; i++) /* Cubature */ for (j = 0; j <= KubY; j++) { *Wert += K4KnotGew (i,j) * f (a + i * Hab, c + j * Hcd); } *Wert *= Hab * Hcd * KubVer * KubVer; /* final multipl.*/ if (Schaetzen && k == 1) Wert1 = *Wert; /* store value */ } if (Schaetzen) /* estimate */ { *FehlerSch = (Wert1 - *Wert) / ((1 << Ordnung) - 1); } return (0); } static REAL K4KnotGew (int i, int j) /*********************************************************************** * Local function for finding node weights * * * * Weight the functional values at the nodes depending on their location* * (edge, interior, ...) in the summed cubature formulas. * ***********************************************************************/ { REAL f; int k; #define Faktor(a,b) ((k == 0 && a > 0 && a < b) ? TWO : ONE) /*********************************************************** * for x-direction : * * 1) node at interval end (k == 0) * * 2) node not at left end (a > 0) nor at right end (a 0) and not at top (a < b) * ***********************************************************/ k = i % KubVer; f = KubArr [KubVer-1] [k] * Faktor (i, KubX); k = j % KubVer; f *= KubArr [KubVer-1] [k] * Faktor (j, KubY); #undef Faktor return (f); } int Kub3NeC3 ( REAL Px, REAL Py, REAL Qx, REAL Qy, REAL Rx, REAL Ry, int n, REAL f (REAL,REAL), REAL* Wert ); /*********************************************************************** * Cubature over triangles using the Newton-Cotes formulas * * * * This function integrates f (x,y) over the triangle PQR using the * * summed 3 point cubature formulas of Newton-Cotes on sub-triangles. * * The sub-triangles are obtained by a regular partition of the edges * * of PQR into n equal parts. * * * * Input parameters: * * REAL Px,Py coordinates of P * * REAL Qx,Ry coordinates of Q * * REAL Rx,Ry coordinates of R * * int n partion number along edges * * REAL f () function * * REAL *Wert value of integral * * * * Return value : * * 0: o.k. * * 1: n improper * * 2: corners P, Q and R are collinear * * * * Author: Uli Eggermann, 8.1.1990 * ***********************************************************************/ #define Kub3NeC3Epsilon .000001 /* for collinearity test */ /*********************************************************************** * Cubature over triangle via summed 3-point Newton-Cotes formula * ***********************************************************************/ int Kub3NeC3 (REAL Px, REAL Py, REAL Qx, REAL Qy, REAL Rx, REAL Ry, int n, REAL f (REAL,REAL), REAL* Wert) { REAL hPQx, hPQy, hPRx, hPRy, Area; int i = 0, j; if (n < 1) return 1; /* n ok ? */ Area = Px * (Qy - Ry) + Qx * (Ry - Py) + Rx * (Py - Qy); /* P, Q and R */ if (FABS (Area) < Kub3NeC3Epsilon ) return 2; /* collinear? */ n *= 2; /* number of halved edges */ hPQx = (Qx - Px) / n; hPQy = (Qy - Py) / n; /* halve the vector PQ */ hPRx = (Rx - Px) / n; hPRy = (Ry - Py) / n; /* halve the vector PR */ #define auswerten (i % 2 || j % 2) /* i or j even */ #define X Px + hPQx*i + hPRx*j /* x-argument */ #define Y Py + hPQy*i + hPRy*j /* y-argument */ #define wieoft ((i==0 || j==0 || i==n-j) ? 1 : 2) /* do not double at the edges */ *Wert = ZERO; /* integral = 0 */ for (j = 0; j < n; j++) /* j moves along PR */ for (i = 0; i <= n-j; i++) /* i moves along PQ */ if (auswerten) /* sum if needed */ *Wert += wieoft * f(X,Y); *Wert += wieoft * f(X,Y); /* sum */ *Wert *= Area / ((REAL)1.5 * n * n); /* last factor */ return 0; #undef Kub3NeC3Epsilon #undef auswerten #undef X #undef Y #undef wieoft } int Kub4RoRi ( REAL a, REAL b, int Nx, REAL c, REAL d, int Ny, int nST, REAL f (REAL,REAL), REAL* Wert, REAL* FehlerSch ) /*********************************************************************** * Cubature over a rectangle using the summed Romberg-Richardson method * * * * Intecrate the function f(x,y) over a rectangle using the summed * * trapezoidal rule for the Romberg sequence of step sizes. Then we * * compute an improved value using Richardson extrapolations. * * * * The Richardson sequence is : * * { h/2, h/4, h/8, h/16, h/32, h/64, h/128, h/256 .. } * * * * Parameters: * * REAL a,b,Nx left, right x-end points * * int Nx # of x-intervals * * REAL c,d,Ny ditto for y * * int Ny # of y-intervals * * int nST number of summed trapezoidal cubatures * * REAL f () function * * REAL *Wert value for integral * * REAL *FehlSch error estimate for Wert * * * * Return value : * * 0: o.k. * * 1: Nx improper * * 2: Ny improper * * 3: nST improper * * 4: rectangle has area zero * * 5: lack of available memory * * * * subroutine used : * * Kub4RoST (lokal) * * * * Author Uli Eggermann, 8.1.1990 * ***********************************************************************/ { int ret; REAL *L; void *vmblock; if (Nx < 1) return (1); if (Ny < 1) return (2); if (nST < 1) return (3); if (a == b || c == d) return (4); vmblock = vminit(); L = (REAL *)vmalloc(vmblock, VEKTOR, nST, 0); if (! vmcomplete(vmblock)) return 5; /* initialize Romberg column : */ /* Cubature using trapezoidal rule */ Kub4RoST (a, b, Nx, c, d, Ny, nST, f, L); /* Compute Richardson columns */ ret = RoRiExtr (L, nST, 2, Wert, FehlerSch) ? 5 : 0; vmfree(vmblock); return (ret); } static void Kub4RoST ( REAL a, REAL b, int Nx, REAL c, REAL d, int Ny, int nST, REAL f (REAL,REAL), REAL* W ) /*********************************************************************** * Initialize Romberg column: Cubature via trapezoidal rule * ***********************************************************************/ { int Nab, Ncd, i, j, k; REAL Hab, Hcd, Faktor; for (j = 0; j < nST; j++) { i = 1 << j; Nab = Nx * i; Hab = (b - a) / Nab; /* i = 2^j */ Ncd = Ny * i; Hcd = (d - c) / Ncd; W[j] = ZERO; for (i=0; i<=Nab; i++) for (k=0; k<=Ncd; k++) { Faktor = ONE; /* at edge: 1 */ if (i>0 && i0 && k 0 && i < Nab ) Faktor = TWO; if ( k > 0 && k < Ncd ) Faktor *= TWO; L[j] += Faktor * f (a + i * Hab, c + k * Hcd); } /* hx hy */ L[j] *= Hab * Hcd * (REAL)0.25; /* Multiply by ----- */ /* 4 */ /* printf ("%"LZP"f \n", L(j)); */ } /* ------------------------------------------------------------------ */ int Kub3RoRi ( REAL Px, REAL Py, REAL Qx, REAL Qy, REAL Rx, REAL Ry, int n, REAL f (REAL,REAL), REAL* Wert, REAL* FehlerSch ) /*********************************************************************** * cubature over triangular regions using the summed 3 point formula * * and Romberg-Richardson extrapolation: * * * * Function f (x,y) is integrated approximately over the triangle PQR * * using the summed 3 point cubature formula on similar subtriangles. * * The number of cubatures with halved subtriangle edge length each is * * defined by n. * * * * Input parameters: * * REAL Px,Py coordinaten of P * * REAL Qx,Ry coordinaten of Q * * REAL Rx,Ry coordinaten of R * * int n number of cubatures * * REAL f (x,y) function to be integrated * * * * Output parameters: * * REAL *Wert approximation of the double integral * * REAL *FehlSch error estimate for `Wert' * * * * Return value: * * 0: o.k. * * 1: n improper * * 2: corners P, Q and R are collinear * * 3: no more memory for auxiliary vector * * * * Subroutines used: * * RoRiExtr Romberg-Richardson extrapolation * * Kub3NeC3 computation of cubature value * * Kub3Nec3n computation and storing of the n cubature values * * * * author Uli Eggermann, 07.05.1990 * ***********************************************************************/ { int i; REAL *W; void *vmblock; if (n < 1) return (1); vmblock = vminit(); W = (REAL *)vmalloc(vmblock, VEKTOR, n, 0); if (! vmcomplete(vmblock)) return 3; /* storing Newton-Cotes values in the first column */ i = Kub3Nec3n (Px,Py, Qx,Qy, Rx,Ry, n, f, W); /* computing the resulting Richardson columns, */ if (i == 0) /* but only if the first column is ok. */ if (n == 1) { *Wert = *W; *FehlerSch = ZERO; i = 0; } else i = RoRiExtr (W, n, 2, Wert, FehlerSch); vmfree(vmblock); /* re-allocation */ return i; } static int Kub3Nec3n ( REAL Px, REAL Py, REAL Qx, REAL Qy, REAL Rx, REAL Ry, int n, REAL f(REAL,REAL), REAL* W ) /*********************************************************************** * computing and storing the n cubature values in given array W[]. * ***********************************************************************/ { int i, ret = 0; if (n < 1) return (1); for (i = 0; i < n; i++) { if ((ret = Kub3NeC3 (Px,Py,Qx,Qy,Rx,Ry, 1< 0; p *= s, k--) for (j = 0; j < k; j++) RiEx[j] = (p * RiEx[j+1] - RiEx[j]) / (p - 1); *Wert = RiEx[0]; *FehlerSch = RiEx[0] - RiEx[1]; vmfree(vmblock); return (0); } /* ------------------------ END kubnec.cpp -------------------------- */