الحاسبة الشاملة — مشروع الهندرة

import numpy as np
from scipy.integrate import quad

def g(u, xp=21.0):
    if u <= 0 or u >= 2*np.pi: return 0.0
    return (u/(2*np.pi))**5 * np.log(2*np.pi/u) * np.exp(-2*np.pi/(xp*u))

def Y(n, u): return np.sin(n*u/2.0)
def D(u): return np.cos(u/2.0)

xp = 21
N1 = quad(lambda u: (g(u,xp)*Y(1,u))**2, 0, 2*np.pi, limit=300)[0]
N2 = quad(lambda u: (g(u,xp)*Y(2,u))**2, 0, 2*np.pi, limit=300)[0]
N3 = quad(lambda u: (g(u,xp)*Y(3,u))**2, 0, 2*np.pi, limit=300)[0]

A12 = quad(lambda u: (g(u,xp)*Y(1,u)/np.sqrt(N1)) * D(u) * (g(u,xp)*Y(2,u)/np.sqrt(N2)), 0, 2*np.pi, limit=300)[0]
A13 = quad(lambda u: (g(u,xp)*Y(1,u)/np.sqrt(N1)) * D(u) * (g(u,xp)*Y(3,u)/np.sqrt(N3)), 0, 2*np.pi, limit=300)[0]

R = (A13/A12)**2
print(f"N1 = {N1:.10e}")
print(f"N2 = {N2:.10e}")
print(f"N3 = {N3:.10e}")
print(f"A12 = {A12:.10e}")
print(f"A13 = {A13:.10e}")
print(f"R   = {R:.10f}")
print(f"Target error = {abs(R - 0.60352):.3e}")

import numpy as np

def g(u, xp=21.0):
    """Envelope function"""
    return (u/(2*np.pi))**5 * np.log(2*np.pi/u) * np.exp(-2*np.pi/(xp*u))

def p(u, xp=21.0):
    """p = g'/g"""
    ell = np.log(2*np.pi/u)
    a = 2*np.pi/xp
    return 5/u - 1/(u*ell) + a/(u**2)

def dp(u, xp=21.0):
    """p' = d/du(g'/g)"""
    ell = np.log(2*np.pi/u)
    a = 2*np.pi/xp
    return -5/u**2 + (ell-1)/(u**2*ell**2) - 2*a/u**3

def g_second_analytic(u, xp=21.0):
    """g'' analytically = g * (p^2 + p')"""
    return g(u,xp) * (p(u,xp)**2 + dp(u,xp))

def g_second_numeric(u, xp=21.0, h=1e-5):
    """g'' numerically via central differences"""
    return (g(u+h,xp) - 2*g(u,xp) + g(u-h,xp)) / (h**2)

u, xp, h = 1.0, 21.0, 1e-5

gu = g(u, xp)
gpp_num = g_second_numeric(u, xp, h)
gpp_ana = g_second_analytic(u, xp)

print(f"u = {u}, xp = {xp}, h = {h}")
print(f"g(u)              = {gu:.10e}")
print(f"g'' (numeric)     = {gpp_num:.10e}")
print(f"g'' (analytic)    = {gpp_ana:.10e}")
print(f"derivative diff   = {abs(gpp_num - gpp_ana):.4e}")
print()
print("Identity check: -g'' + (g''/g)*g = 0, where g != 0.")
print("This is an algebraic identity, not a spectral verification.")