' Numeric Integration — rectangle, trapezoid, Simpson rules
#COMPILE EXE
#DIM ALL

FUNCTION F1(x AS DOUBLE) AS DOUBLE
    FUNCTION = x * x
END FUNCTION

FUNCTION F2(x AS DOUBLE) AS DOUBLE
    FUNCTION = SIN(x)
END FUNCTION

FUNCTION F3(x AS DOUBLE) AS DOUBLE
    FUNCTION = EXP(-x * x)
END FUNCTION

FUNCTION RectangleRule(func AS LONG, a AS DOUBLE, b AS DOUBLE, n AS LONG) AS DOUBLE
    LOCAL h AS DOUBLE
    LOCAL sum AS DOUBLE
    LOCAL i AS LONG
    LOCAL x AS DOUBLE
    h = (b - a) / n
    sum = 0.0
    FOR i = 0 TO n - 1
        x = a + (i + 0.5) * h
        IF func = 1 THEN sum = sum + F1(x)
        IF func = 2 THEN sum = sum + F2(x)
        IF func = 3 THEN sum = sum + F3(x)
    NEXT i
    FUNCTION = sum * h
END FUNCTION

FUNCTION TrapezoidRule(func AS LONG, a AS DOUBLE, b AS DOUBLE, n AS LONG) AS DOUBLE
    LOCAL h AS DOUBLE
    LOCAL sum AS DOUBLE
    LOCAL i AS LONG
    LOCAL x AS DOUBLE
    h = (b - a) / n
    IF func = 1 THEN sum = (F1(a) + F1(b)) / 2.0
    IF func = 2 THEN sum = (F2(a) + F2(b)) / 2.0
    IF func = 3 THEN sum = (F3(a) + F3(b)) / 2.0
    FOR i = 1 TO n - 1
        x = a + i * h
        IF func = 1 THEN sum = sum + F1(x)
        IF func = 2 THEN sum = sum + F2(x)
        IF func = 3 THEN sum = sum + F3(x)
    NEXT i
    FUNCTION = sum * h
END FUNCTION

FUNCTION SimpsonRule(func AS LONG, a AS DOUBLE, b AS DOUBLE, n AS LONG) AS DOUBLE
    LOCAL h AS DOUBLE
    LOCAL sum AS DOUBLE
    LOCAL i AS LONG
    LOCAL x AS DOUBLE
    IF n MOD 2 = 1 THEN n = n + 1
    h = (b - a) / n
    IF func = 1 THEN sum = F1(a) + F1(b)
    IF func = 2 THEN sum = F2(a) + F2(b)
    IF func = 3 THEN sum = F3(a) + F3(b)
    FOR i = 1 TO n - 1
        x = a + i * h
        IF i MOD 2 = 0 THEN
            IF func = 1 THEN sum = sum + 2.0 * F1(x)
            IF func = 2 THEN sum = sum + 2.0 * F2(x)
            IF func = 3 THEN sum = sum + 2.0 * F3(x)
        ELSE
            IF func = 1 THEN sum = sum + 4.0 * F1(x)
            IF func = 2 THEN sum = sum + 4.0 * F2(x)
            IF func = 3 THEN sum = sum + 4.0 * F3(x)
        END IF
    NEXT i
    FUNCTION = sum * h / 3.0
END FUNCTION

SUB CompareMethods()
    LOCAL a AS DOUBLE
    LOCAL b AS DOUBLE
    LOCAL n AS LONG
    LOCAL rect AS DOUBLE
    LOCAL trap AS DOUBLE
    LOCAL simp AS DOUBLE
    a = 0.0
    b = 1.0
    n = 100
    PRINT "=== COMPARISON: Integral of x^2 from 0 to 1 ==="
    PRINT "Exact value: 0.333333..."
    PRINT "N=" & STR$(n)
    rect = RectangleRule(1, a, b, n)
    trap = TrapezoidRule(1, a, b, n)
    simp = SimpsonRule(1, a, b, n)
    PRINT "Rectangle: " & FORMAT$(rect, "##.########") & " (error: " & FORMAT$(rect - 1.0/3.0, "##.########") & ")"
    PRINT "Trapezoid: " & FORMAT$(trap, "##.########") & " (error: " & FORMAT$(trap - 1.0/3.0, "##.########") & ")"
    PRINT "Simpson:   " & FORMAT$(simp, "##.########") & " (error: " & FORMAT$(simp - 1.0/3.0, "##.########") & ")"
    PRINT
END SUB

SUB ConvergenceStudy()
    LOCAL n AS LONG
    LOCAL trap AS DOUBLE
    LOCAL exact AS DOUBLE
    a = 0.0
    b = 1.0
    exact = 1.0 / 3.0
    PRINT "=== CONVERGENCE: Trapezoid Rule ==="
    PRINT "Integral of x^2 from 0 to 1"
    PRINT "N      Result            Error"
    PRINT STRING$(50, 61)
    n = 10
    DO WHILE n <= 1000
        trap = TrapezoidRule(1, a, b, n)
        PRINT FORMAT$(n, "#####") & "  " & FORMAT$(trap, "##.##########") & "      " & FORMAT$(ABS(trap - exact), "##.##########")
        n = n * 2
    LOOP
    PRINT
END SUB

SUB IntegrateSine()
    LOCAL n AS LONG
    LOCAL simp AS DOUBLE
    a = 0.0
    b = 3.14159265358979
    PRINT "=== INTEGRAL OF SIN(X) FROM 0 TO PI ==="
    PRINT "Exact: 2.0"
    PRINT "N      Simpson"
    PRINT STRING$(30, 61)
    n = 10
    DO WHILE n <= 160
        simp = SimpsonRule(2, a, b, n)
        PRINT FORMAT$(n, "#####") & "  " & FORMAT$(simp, "##.##########")
        n = n * 2
    LOOP
    PRINT
END SUB

SUB GaussianIntegral()
    LOCAL n AS LONG
    LOCAL rect AS DOUBLE
    a = -2.0
    b = 2.0
    PRINT "=== INTEGRAL OF EXP(-X^2) FROM -2 TO 2 ==="
    PRINT "Approximate exact: SQRT(PI) * erf(2) ~ 1.76416"
    PRINT "N      Rectangle"
    n = 100
    rect = RectangleRule(3, a, b, n)
    PRINT FORMAT$(n, "#####") & "  " & FORMAT$(rect, "##.##########")
    n = 1000
    rect = RectangleRule(3, a, b, n)
    PRINT FORMAT$(n, "#####") & "  " & FORMAT$(rect, "##.##########")
    PRINT
END SUB

FUNCTION PBMAIN() AS LONG
    PRINT "============================="
    PRINT " Numeric Integration Demo"
    PRINT "============================="
    PRINT
    CALL CompareMethods()
    CALL ConvergenceStudy()
    CALL IntegrateSine()
    CALL GaussianIntegral()
    PRINT "Integration demos complete."
    FUNCTION = 0
END FUNCTION
