DECLARE SUB IntegralSimpson (kind%, a#, b#, n%, res#) DECLARE FUNCTION F# (kind%, t#) DECLARE FUNCTION sinintegral# (x#) DECLARE FUNCTION cosintegral# (x#) '************************************************************** '* Program to demonstrate the use of functions sinintegral(x) * '* or cosintegral(x) using the Simpson's method * '* ---------------------------------------------------------- * '* * '* SAMPLE RUNS: * '* * '* Kind (1 or 2)= 1 * '* * '* Test of Function sinintegral(x) by Simpson's method * '* * '* Input value for x variable: * '* * '* X = 2 * '* * '* Result: 1.605412976802747 * '* * '* X = 5 * '* * '* Result: 1.549931244944637 * '* * '* X = 80 * '* * '* Result: 1.572330886912482 * '* * '* * '* Kind (1 or 2) = 2 * '* * '* Test of Function cosintegral(x) by Simpson's method * '* * '* Input value for x variable: * '* * '* X = 2 * '* * '* Result: .4229808287748087 * '* * '* X = 5 * '* * '* Result: -.1900297496567194 * '* * '* X = 80 * '* * '* Result: -1.240250115512912D-02 * '* * '* Basic Version By J-P Moreau. * '* (www.jpmoreau.fr) * '************************************************************** 'Program Test_Sinintegral DEFDBL A-H, O-Z DEFINT I-N CLS PRINT INPUT " Kind (1 or 2) = ", kind PRINT IF kind = 1 THEN PRINT " Test of Function sinintegral(x) by Simpson's method" ELSE PRINT " Test of Function cosintegral(x) by Simpson's method" END IF PRINT PRINT " Input value for x variable:" PRINT INPUT " X = ", x PRINT IF kind = 2 AND x <= 0# THEN PRINT " WARNING! X must be positive." ELSE IF kind = 1 THEN PRINT " Result: "; sinintegral(x) ELSE PRINT " Result: "; cosintegral(x) END IF END IF PRINT END 'of main program '******************************************************* '* Calcualte Ci(x)= gamma + Ln(x) + I(x) with I(x) = * '* Integral of (cos(u)-1)/u from u=0 to u=x (x > 0). * '* Gamma = 0.57721566... (Euler's constant). * '******************************************************* FUNCTION cosintegral (x) 'x0 begin x value 'x end x value 'nstep number of integration steps 'res result of integral gamma = .5772156649015329# 'Euler's constant x0 = 0# nstep = INT(200 * ABS(x)) 'this should ensure about 14 exact digits IntegralSimpson 2, x0, x, nstep, res 'kind=2 cosintegral = gamma + LOG(x) + res END FUNCTION 'Given function to integrate FUNCTION F (kind, t) IF kind = 1 THEN IF ABS(t) < .0000000001# THEN F = 1# ELSE F = SIN(t) / t 'for sinintegral END IF ELSE IF ABS(t) < .0000000001# THEN F = 0# ELSE F = (COS(t) - 1#) / t 'for cosintegral END IF END IF END FUNCTION '******************************************************* '* Integral of a function F(X) by Simpson's method * '* --------------------------------------------------- * '* INPUTS: * '* kind =1 for sinintegral * '* =2 for cosintegral * '* a begin value of x variable * '* b end value of x variable * '* n number of integration steps * '* * '* OUTPUT: * '* res the integral of FUNC(X) from a to b * '* * '******************************************************* SUB IntegralSimpson (kind, a, b, n, res) xstep = (b - a) / 2 / n r = F(kind, a) res = (r + F(kind, b)) / 2# FOR i = 1 TO 2 * n - 1 r = F(kind, a + i * xstep) IF (i MOD 2) <> 0 THEN res = res + r + r ELSE res = res + r END IF NEXT i res = res * xstep * 2# / 3# END SUB '******************************************************* '* Si(x) = Integral of sin(u)/u from u=0 to u=x * '******************************************************* FUNCTION sinintegral (x) 'x0 begin x value 'x end x value 'nstep number of integration steps 'res result of integral x0 = 0# nstep = INT(200 * ABS(x)) 'this should ensure about 14 exact digits IntegralSimpson 1, x0, x, nstep, res 'kind=1 sinintegral = res END FUNCTION