'*********************************************************** '* Test program for cubature over triangles using 3-point * '* Newton-Cotes and Romberg-Richardson extrapolation * '* ------------------------------------------------------- * '* SAMPLE RUN: * '* (Integrate function f(x,y)=xy in triangle defined by * '* 3 points (0, 0), (10, 0) and (0, 10) ). * '* * '* ------------------------------------------------------- * '* # Newton- Error Value Error * '* Cotes Code Integral Estimation * '* ------------------------------------------------------- * '* 1 0 4166.666667 0 * '* 2 0 777.777778 -185.7638888888889 * '* 3 0 507.936508 -2.480158730158735 * '* 4 0 451.557907 -.2167842445620067 * '* 5 0 432.920438 -1.887945480257258D-02 * '* 6 0 424.647389 -2.020102657809275D-03 * '* 7 0 420.638618 -2.446763548391573D-04 * '* 8 0 418.650342 -3.033867784552058D-05 * '* 9 0 417.658217 -3.784657792493817D-06 * '* 10 0 417.162406 -4.72842373255844D-07 * '* 11 0 416.914532 -5.90977720094088D-08 * '* 12 0 416.790599 -7.386972811218584D-09 * '* ------------------------------------------------------- * '* Reference: * '* "Numerical Algorithms with C, By Gisela Engeln-Muellges * '* and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * '* * '* Basic Release 1.0 By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '*********************************************************** 'Program TK3NEC DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 DIM W(25) AS DOUBLE DIM RiEx(25) AS DOUBLE F$ = "#####.######" F1$ = " ##" CLS PRINT PRINT "-------------------------------------------------------" PRINT " # Newton- Error Value Error " PRINT " Cotes Code Integral Estimation " PRINT "-------------------------------------------------------" FOR j = 1 TO 12 Px = 0#: Py = 0#: Qx = 10#: Qy = 0#: Rx = 0#: Ry = 10#: n = j GOSUB 1000 'call Kub3RoRi(Px,Py,Qx,Qy,Rx,Ry, n, Wert, FehlerSch, icode) PRINT USING F1$; j; PRINT " "; icode; " "; PRINT USING F$; Wert; PRINT " "; FehlerSch NEXT j PRINT "-------------------------------------------------------" END 'of main program 'User defined function 500 'func(x,y) func = x * y RETURN 1000 'Subroutine Kub3RoRi(Px, Py, Qx, Qy, Rx, Ry, n, Wert, FehlerSch,icode) '************************************************************************ '* 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: * '* Double Px,Py coordinaten of P * '* Double Qx,Ry coordinaten of Q * '* Double Rx,Ry coordinaten of R * '* integer n number of cubatures * '* * '* Output parameters: * '* Double Wert approximation of the double integral * '* Double FehlSch error estimate for Wert * '* * '* Return value (icode): * '* 0: o.k. * '* 1: n improper * '* 2: corners P, Q and R are collinear * '* 3: no more memory for auxiliary vector * '* (not used here). * '* * '* 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 (C version) * '************************************************************************ IF n < 1 THEN icode = 1 RETURN END IF 'storing Newton-Cotes values in the first column GOSUB 3000 'call Kub3Nec3n (Px,Py, Qx,Qy, Rx,Ry, n, W,icode) 'computing the resulting Richardson columns, IF icode = 0 THEN 'but only if the first column is ok. IF n = 1 THEN Wert = W(0) FehlerSch = 0# icode = 0 ELSE Ordnung = 2: n2 = n GOSUB 4000 'call RoRiExtr (W, n2, ordnung, Wert, FehlerSch,icode) END IF END IF RETURN 2000 'Subroutine Kub3NeC3(Px, Py, Qx, Qy, Rx, Ry, n1, Wert, ic) '************************************************************************ '* 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: * '* Double Px,Py coordinates of P * '* Double Qx,Ry coordinates of Q * '* Double Rx,Ry coordinates of R * '* integer n1 partion number along edges * '* Double Wert value of integral * '* * '* Return value ic: * '* 0: o.k. * '* 1: n1 improper * '* 2: corners P, Q and R are collinear * '* * '* Author: Uli Eggermann, 8.1.1990 (C version) * '************************************************************************ '************************************************************************ '* Cubature over triangle via summed 3-point Newton-Cotes formula * '************************************************************************ DIM wieoft AS INTEGER IF n1 < 1 THEN ic = 1 RETURN 'n1 ok ? END IF Area = Px * (Qy - Ry) + Qx * (Ry - Py) + Rx * (Py - Qy) 'P, Q and R IF ABS(Area) < .000001 THEN 'collinear? ic = 2 RETURN END IF n1 = n1 * 2 'number of halved edges hPQx = (Qx - Px) / n1 hPQy = (Qy - Py) / n1 'halve the vector PQ hPRx = (Rx - Px) / n1 hPRy = (Ry - Py) / n1 'halve the vector PR Wert = 0# 'integral = 0 FOR jj = 0 TO n1 - 1 'jj moves along PR FOR ii = 0 TO n1 - jj 'ii moves along PQ IF (ii MOD 2) <> 0 OR (jj MOD 2) <> 0 THEN 'sum is needed x = Px + hPQx * ii + hPRx * jj y = Py + hPQy * ii + hPRy * jj GOSUB 500 'calculate func(x,y) IF ii = 0 OR jj = 0 OR ii = n - jj THEN wieoft = 1 ELSE wieoft = 2 END IF Wert = Wert + wieoft * func END IF NEXT ii NEXT jj x = Px + hPQx * ii + hPRx * jj y = Py + hPQy * ii + hPRy * jj GOSUB 500 IF ii = 0 OR jj = 0 OR ii = n - jj THEN wieoft = 1 ELSE wieoft = 2 END IF Wert = Wert + wieoft * func 'sum Wert = Wert * Area / (1.5# * n1 * n1) 'last factor ic = 0 RETURN 3000 'Subroutine Kub3Nec3n(Px, Py, Qx, Qy, Rx, Ry, n, W,icode) '************************************************************************ '* computing and storing the n cubature values in given array W. * '************************************************************************ IF n < 1 THEN icode = 1 RETURN END IF FOR i = 0 TO n - 1 n1 = 2 ^ i GOSUB 2000 'call Kub3NeC3(Px,Py,Qx,Qy,Rx,Ry, 2^i, W(i),ic) W(i) = Wert IF ic <> 0 THEN GOTO 10 NEXT i 10 icode = ic RETURN 4000 'Subroutine RoRiExtr(W, n2, Ordnung, Wert, FehlerSch,icode) '************************************************************************ '* Richardson extrapolation for a given Romberg sequence * '* * '* Parameter: * '* Double W(n) given Romberg sequence * '* integer n2 length of W * '* integer Ordnung Order of method * '* Double Wert final value of extrapolation * '* Double FehlerSch error estimate of Wert * '* * '* Return value icode: * '* 0: o.k. * '* 1: lack of available memory (not used here) * '* 2: n2 too small * '* * '* Author Uli Eggermann, 3.31.1996 * '************************************************************************ IF n2 < 2 THEN icode = 2 RETURN END IF FOR k = 0 TO n2 - 1 'copy RiEx(k) = W(k) NEXT k s = 2 ^ Ordnung p = s FOR k = n2 - 1 TO 1 STEP -1 p = p * s FOR jj = 0 TO k - 1 RiEx(jj) = (p * RiEx(jj + 1) - RiEx(jj)) / (p - 1) NEXT jj NEXT k Wert = RiEx(0) FehlerSch = RiEx(0) - RiEx(1) icode = 0 RETURN 'end of file tk3nec.bas