zerr.py

1import mpmath as mp2import sympy as sp34DPS = 505SAMPLES = 5067def symbolic_residuals():8    a, u = sp.symbols("a u", positive=True)9    r = 1 - a**2 + a**2 * u**210    cube = (-a * u + sp.sqrt(r)) ** 311    t1 = -(a**3) * u**312    t2 = 3 * a**2 * u**2 * sp.sqrt(r)13    t3 = -3 * a * u * r14    t4 = r ** sp.Rational(3, 2)15    even = sp.sqrt(r) * (1 - a**2 + 4 * a**2 * u**2)16    prim = u * r ** sp.Rational(3, 2)17    res = {18        "cube minus its four terms": sp.simplify(sp.expand(cube - (t1 + t2 + t3 + t4))),19        "even part minus sqrt(R)(1-a^2+4a^2u^2)": sp.simplify(even - (t2 + t4)),20        "d/du[u R^(3/2)] minus even part": sp.simplify(sp.diff(prim, u) - even),21    }22    boundary = sp.simplify(prim.subs(u, 1) - prim.subs(u, -1))23    closed = sp.simplify(sp.integrate(cube, (u, -1, 1)))24    return res, boundary, closed2526def inner_integral(a):27    return mp.quad(lambda u: (-a * u + mp.sqrt(1 - a**2 + a**2 * u**2)) ** 3, [-1, 0, 1])2829def chord(a, phi):30    return -a * mp.cos(phi) + mp.sqrt(1 - a**2 * mp.sin(phi) ** 2)3132def i_diam_direct():33    half = mp.pi / 234    return mp.quad(35        lambda a, phi: chord(a, phi) ** 3 * mp.sin(phi),36        [-1, 0, 1],37        [0, half, mp.pi],38    )3940def main():41    mp.mp.dps = DPS42    print("ZERR DECOMPOSITION")43    res, boundary, closed = symbolic_residuals()44    for name, value in res.items():45        print(f"  residual {name} = {value}")46    print(f"  boundary [u R^(3/2)] from -1 to 1 = {boundary}")47    print(f"  closed form of the inner integral = {closed}")4849    worst = mp.mpf(0)50    for i in range(SAMPLES):51        a = mp.mpf(i + 1) / (SAMPLES + 1)52        worst = max(worst, abs(inner_integral(a) - 2))53    print(f"  inner integral at {SAMPLES} values of a, working digits {DPS}")54    print(f"  max deviation from 2 = {mp.nstr(worst, 5)}")5556    i_diam = mp.quad(inner_integral, [-1, 0, 1])57    print(f"  I_diam by the constant inner integral = {mp.nstr(i_diam, 30)}")5859    mp.mp.dps = 2560    direct = i_diam_direct()61    print(f"  I_diam by direct chord-cube quadrature = {mp.nstr(direct, 20)}")6263    mp.mp.dps = DPS64    area = mp.pi / 265    prob = mp.mpf(4) / (3 * area**2)66    exact = 16 / (3 * mp.pi**2)67    print(f"  Area(H) = pi/2, Area(H)^2 = pi^2/4 = {mp.nstr(area**2, 20)}")68    print(f"  P(chord crosses the diameter) = I_diam/(3 Area(H)^2) = {mp.nstr(prob, 20)}")69    print(f"  16/(3 Pi^2) = {mp.nstr(exact, 20)}")70    print(f"  16/(3 Pi^2) to 16 places = {mp.nstr(exact, 16)}")71    print(f"  1 - 16/(3 Pi^2) to 16 places = {mp.nstr(1 - exact, 16)}")72    print(f"  |P - 16/(3 Pi^2)| = {mp.nstr(abs(prob - exact), 5)}")7374main()