কম্পোজিট রুল ও অ্যাডাপটিভ কোয়াড্রেচার

import math

def f(x):
    return 1.0 / (1.0 + 1000*(x-5.0)**2)

a, b = 0.0, 10.0
sa = math.sqrt(1000)
true_val = (1.0/sa) * (math.atan(sa*(b-5.0)) - math.atan(sa*(a-5.0)))
print(f"true integral = {true_val:.10f}")
print()

def simpson(f, a, b, n):
    h = (b-a)/n
    total = f(a) + f(b)
    for i in range(1, n):
        x = a + i*h
        total += (4 if i % 2 != 0 else 2)*f(x)
    return total*h/3

print("স্থির/ইউনিফর্ম গ্রিডে কম্পোজিট সিম্পসনস:")
for n in [16, 64, 256, 1024]:
    approx = simpson(f, a, b, n)
    err = abs(approx-true_val)
    print(f"  n={n:>5} (evals={n+1:>5})  approx={approx:.10f}  error={err:.3e}")

import math

def f(x):
    return 1.0 / (1.0 + 1000*(x-5.0)**2)

a, b = 0.0, 10.0
sa = math.sqrt(1000)
true_val = (1.0/sa) * (math.atan(sa*(b-5.0)) - math.atan(sa*(a-5.0)))

def simpson_rule(fa, fm, fb, a, b):
    return (b - a) / 6 * (fa + 4*fm + fb)

def adaptive_simpson(f, a, b, tol, fa=None, fm=None, fb=None, depth=0, max_depth=50):
    if fa is None: fa = f(a)
    if fm is None: fm = f((a+b)/2)
    if fb is None: fb = f(b)
    S = simpson_rule(fa, fm, fb, a, b)

    m = (a+b)/2
    lm, rm = (a+m)/2, (m+b)/2
    flm, frm = f(lm), f(rm)
    S_left  = simpson_rule(fa, flm, fm, a, m)
    S_right = simpson_rule(fm, frm, fb, m, b)
    S2 = S_left + S_right

    if depth >= max_depth or abs(S2 - S) < 15*tol:
        return S2 + (S2 - S)/15   # রিচার্ডসন-স্টাইল সংশোধন, L22 দেখুন
    return (adaptive_simpson(f, a, m, tol/2, fa, flm, fm, depth+1, max_depth) +
            adaptive_simpson(f, m, b, tol/2, fm, frm, fb, depth+1, max_depth))

evals_count = [0]
def f_counted(x):
    evals_count[0] += 1
    return f(x)

result = adaptive_simpson(f_counted, a, b, 1e-8)
err = abs(result - true_val)
print(f"true integral        = {true_val:.10f}")
print(f"adaptive Simpson     = {result:.10f}")
print(f"error                = {err:.3e}")
print(f"function evaluations = {evals_count[0]}")