'****************************************************************** '* Test program for cubature over rectangles using Gauss * '* * '* -------------------------------------------------------------- * '* SAMPLE RUN: * '* (Integrate function EXP(-(x*x + y*y)) in rectangle defined by * '* x1=0, x2=10, y1=0, y2=10). * '* * '* -------------------------------------------------------------- * '* # Error Value Error Function * '* Method code Integral Estimation calls * '* -------------------------------------------------------------- * '* 0 0 .2746078818715234 16 * '* 1 0 .9434560943928603 64 * '* 2 0 .7739036744095703 144 * '* 3 0 .7853778748173256 256 * '* 4 0 .7854721566819686 400 * '* 5 0 .7853906291687811 576 * '* 6 0 .7853986174112774 784 * '* 7 0 .785398144613683 1024 * '* 0 0 .7797341623610217 .1683754268298328 80 * '* 1 0 .7867640862777988 -5.223066937168715D-02 320 * '* 2 0 .7852707578415222 3.789027810650639D-03 720 * '* 3 0 .7854042225690866 8.78258392034148D-06 1280 * '* 4 0 .7853979886123033 -2.472268988843821D-05 2000 * '* 5 0 .7853981667851698 2.512538796239916D-06 2880 * '* 6 0 .7853981633502606 -1.513536722619335D-07 3920 * '* 7 0 .7853981633979421 6.261419716047101D-09 5120 * '* -------------------------------------------------------------- * '* 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 TK4GAU DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 DIM Verfahren AS INTEGER Nx = 4: Ny = 4 a = 0#: c = 0# b = 10#: d = 10# CLS PRINT " --------------------------------------------------------------" PRINT " # Error Value Error Function " PRINT " Method code Integral Estimation calls " PRINT " --------------------------------------------------------------" FOR mSCH = 0 TO 1 FOR Verfahren = 0 TO 7 ncalls = 0 GOSUB 1000 'call Kub4GauE(a, b, Nx, c, d, Ny, Verfahren, Wert, mSCH, F, ncalls,icode) IF mSCH <> 0 THEN PRINT " "; Verfahren; " "; icode; " "; Wert; F; ELSE PRINT " "; Verfahren; " "; icode; " "; Wert; END IF PRINT TAB(57); USING "####"; ncalls NEXT Verfahren NEXT mSCH PRINT " --------------------------------------------------------------" END 'of main program 'user defined function(x,y) to integrate 500 func = EXP(-(x * x + y * y)) RETURN 1000 'Kub4GauE(a,b,Nx,c,d,Ny,Verf,Wert,Schaetzen,FehlerSch,ncalls,icode) DIM Schaetzen AS INTEGER DIM Verf AS INTEGER Verf = Verfahren: Schaetzen = mSCH '************************************************************************ '* Cubature over rectangles using Newton-Cotes formulas. * '* * '* Integrate the function f(x,y) using the summed cubature formula * '* of Gauss for the rectangle (a,b; c,d). * '* The edges of the sub-rectangles have the lengths: (b-a) / Nx and * '* (d-c) / Ny . * '* * '* Input parameters: * '* DOUBLE a, b, Nx left, right x-end points, number of intervals * '* DOUBLE c, d, Ny left, right y-end points, number of intervals * '* integer Verf order of method (0 <= Verf <= 7) * '* DOUBLE Wert value of integral * '* integer Schaetzen if nonzero : find estimate * '* DOUBLE F error estimate, obtained if desired by a second * '* cubature pass using half of the step size * '* integer ncalls number of function calls * '* * '* Return value icode: * '* 0: o.k. * '* 1: Nx improper * '* 2: Ny improper * '* 3: order incorrect * '* 4: Integration interval has length zero * '* * '* Author Uli Eggermann, 03.31.1996 (C version) * '************************************************************************ DIM u AS INTEGER DIM v AS INTEGER DIM K0(2) AS DOUBLE DIM K1(2, 2) AS DOUBLE DIM K2(3, 2) AS DOUBLE DIM K3(4, 2) AS DOUBLE DIM K4(5, 2) AS DOUBLE DIM K5(6, 2) AS DOUBLE DIM K6(7, 2) AS DOUBLE DIM K7(8, 2) AS DOUBLE ' Initialize K0 to K7 constants K0(0) = 0#: K0(1) = 2# K1(0, 0) = -.577350269189626#: K1(0, 1) = 1! K1(1, 0) = .577350269189626#: K1(1, 1) = 1! K2(0, 0) = -.774596669241483#: K2(0, 1) = .5555555555555556# K2(1, 0) = 0#: K2(1, 1) = .8888888888888888# K2(2, 0) = .774596669241483#: K2(2, 1) = .5555555555555556# K3(0, 0) = -.861136311594053#: K3(0, 1) = .347854845137454# K3(1, 0) = -.339981043584856#: K3(1, 1) = .652145154862546# K3(2, 0) = .339981043584856#: K3(2, 1) = .652145154862546# K3(3, 0) = .861136311594053#: K3(3, 1) = .347854845137454# K4(0, 0) = -.906179845938664#: K4(0, 1) = .236926885056189# K4(1, 0) = -.538469310105683#: K4(1, 1) = .478628670499366# K4(2, 0) = 0#: K4(2, 1) = .5688888888888889# K4(3, 0) = .538469310105683#: K4(3, 1) = .478628670499366# K4(4, 0) = .906179845938664#: K4(4, 1) = .236926885056189# K5(0, 0) = -.9324695142031521#: K5(0, 1) = .17132449237917# K5(1, 0) = -.661209386466265#: K5(1, 1) = .360761573048139# K5(2, 0) = -.238619186083197#: K5(2, 1) = .467913934572691# K5(3, 0) = .238619186083197#: K5(3, 1) = .467913934572691# K5(4, 0) = .661209386466265#: K5(4, 1) = .360761573048139# K5(5, 0) = .9324695142031521#: K5(5, 1) = .17132449237917# K6(0, 0) = -.949107912342759#: K6(0, 1) = .12948496616887# K6(1, 0) = -.741531185599394#: K6(1, 1) = .279705391489277# K6(2, 0) = -.405845151377397#: K6(2, 1) = .381830050505119# K6(3, 0) = 0#: K6(3, 1) = .417959183673469# K6(4, 0) = .405845151377397#: K6(4, 1) = .381830050505119# K6(5, 0) = .741531185599394#: K6(5, 1) = .279705391489277# K6(6, 0) = .949107912342759#: K6(6, 1) = .12948496616887# K7(0, 0) = -.960289856497536#: K7(0, 1) = .101228536290376# K7(1, 0) = -.7966664774136269#: K7(1, 1) = .222381034453374# K7(2, 0) = -.525532409916329#: K7(2, 1) = .313706645877887# K7(3, 0) = -.18343464249565#: K7(3, 1) = .362683783378362# K7(4, 0) = .18343464249565#: K7(4, 1) = .362683783378362# K7(5, 0) = .525532409916329#: K7(5, 1) = .313706645877887# K7(6, 0) = .7966664774136269#: K7(6, 1) = .222381034453374# K7(7, 0) = .960289856497536#: K7(7, 1) = .101228536290376# Wert1 = 0# IF Nx < 1 THEN icode = 1 RETURN END IF IF Ny < 1 THEN icode = 2 RETURN END IF IF Verf < 0 OR Verf > 7 THEN icode = 3 RETURN END IF IF a = b OR c = d THEN icode = 4 RETURN END IF IF Schaetzen <> 0 THEN kend = 2 ELSE kend = 1 END IF FOR k = 1 TO kend Nab = k * Nx 'number of intervals Ncd = k * Ny Hab = .5# * (b - a) / Nab 'half of x-step size Hcd = .5# * (d - c) / Ncd 'half of y-step size Wert = 0# 'initialize FOR i = 0 TO Nab - 1 FOR j = 0 TO Ncd - 1 FOR u = 0 TO Verf FOR v = 0 TO Verf IF Verf = 0 THEN w = Hab * Hcd * K0(1) * K0(1) x = a + Hab * (K0(0) + 2 * i + 1) y = c + Hcd * (K0(0) + 2 * j + 1) ELSEIF Verf = 1 THEN w = Hab * Hcd * K1(u, 1) * K1(v, 1) x = a + Hab * (K1(u, 0) + 2 * i + 1) y = c + Hcd * (K1(v, 0) + 2 * j + 1) ELSEIF Verf = 2 THEN w = Hab * Hcd * K2(u, 1) * K2(v, 1) x = a + Hab * (K2(u, 0) + 2 * i + 1) y = c + Hcd * (K2(v, 0) + 2 * j + 1) ELSEIF Verf = 3 THEN w = Hab * Hcd * K3(u, 1) * K3(v, 1) x = a + Hab * (K3(u, 0) + 2 * i + 1) y = c + Hcd * (K3(v, 0) + 2 * j + 1) ELSEIF Verf = 4 THEN w = Hab * Hcd * K4(u, 1) * K4(v, 1) x = a + Hab * (K4(u, 0) + 2 * i + 1) y = c + Hcd * (K4(v, 0) + 2 * j + 1) ELSEIF Verf = 5 THEN w = Hab * Hcd * K5(u, 1) * K5(v, 1) x = a + Hab * (K5(u, 0) + 2 * i + 1) y = c + Hcd * (K5(v, 0) + 2 * j + 1) ELSEIF Verf = 6 THEN w = Hab * Hcd * K6(u, 1) * K6(v, 1) x = a + Hab * (K6(u, 0) + 2 * i + 1) y = c + Hcd * (K6(v, 0) + 2 * j + 1) ELSEIF Verf = 7 THEN w = Hab * Hcd * K7(u, 1) * K7(v, 1) x = a + Hab * (K7(u, 0) + 2 * i + 1) y = c + Hcd * (K7(v, 0) + 2 * j + 1) END IF GOSUB 500: ncalls = ncalls + 1 Wert = Wert + w * func NEXT v NEXT u NEXT j NEXT i IF Schaetzen <> 0 AND k = 1 THEN Wert1 = Wert 'store value NEXT k IF Schaetzen <> 0 THEN 'estimate F = (Wert - Wert1) / 3# END IF Kub4GauE = 0 RETURN 'end of file tk4gau.bas