ch204 — Numerical Limits

import numpy as np
import matplotlib.pyplot as plt

def numerical_limit(f, a, direction='both', steps=20, start_h=1.0):
    """
    Estimate lim_{x->a} f(x) by evaluating f at points approaching a.
    
    Returns: (h_values, f_values) for the approach.
    """
    h_values = start_h / (2.0 ** np.arange(steps))
    results = {}
    
    if direction in ('right', 'both'):
        results['right'] = [(h, f(a + h)) for h in h_values]
    if direction in ('left', 'both'):
        results['left']  = [(h, f(a - h)) for h in h_values]
    
    return results


# Example 1: lim_{x->0} sin(x)/x
f1 = lambda x: np.sin(x) / x
r1 = numerical_limit(f1, 0.0)

print('lim_{x->0} sin(x)/x')
print(f'{"h":>12}  {"from right":>15}  {"from left":>15}')
print('-' * 48)
for (h_r, v_r), (h_l, v_l) in zip(r1['right'], r1['left']):
    print(f'{h_r:>12.2e}  {v_r:>15.10f}  {v_l:>15.10f}')

# Example 2: lim_{x->1} (x^2 - 1)/(x - 1)
# Algebraically this simplifies to x+1, so the limit should be 2
f2 = lambda x: (x**2 - 1) / (x - 1)
r2 = numerical_limit(f2, 1.0)

print('lim_{x->1} (x^2 - 1)/(x - 1)   [expect: 2.0]')
print(f'{"h":>12}  {"from right":>15}  {"from left":>15}')
print('-' * 48)
for (h_r, v_r), (h_l, v_l) in zip(r2['right'][:12], r2['left'][:12]):
    print(f'{h_r:>12.2e}  {v_r:>15.10f}  {v_l:>15.10f}')

# Watch the breakdown: compute (x^2-1)/(x-1) as x -> 1
h_vals = np.array([10**(-k) for k in range(1, 17)])
a = 1.0
f3 = lambda x: (x**2 - 1) / (x - 1)

estimates = []
for h in h_vals:
    try:
        val = f3(a + h)
    except ZeroDivisionError:
        val = np.nan
    estimates.append(val)

errors = [abs(e - 2.0) if not np.isnan(e) else np.nan for e in estimates]

print(f'{"h":>12}  {"estimate":>16}  {"error":>14}')
print('-' * 48)
for h, est, err in zip(h_vals, estimates, errors):
    err_str = f'{err:.2e}' if err is not None and not np.isnan(err) else 'NaN'
    est_str = f'{est:.10f}' if not np.isnan(est) else 'NaN'
    print(f'{h:>12.2e}  {est_str:>16}  {err_str:>14}')

print()
print('For h < 1e-8, floating-point error dominates. The estimate degrades.')
print('This is the same issue seen in ch038 — Precision and Floating Point Errors.')

# Visualize the breakdown
fig, axes = plt.subplots(1, 2, figsize=(13, 5))

valid_h    = [h for h, e in zip(h_vals, errors) if e is not None and not np.isnan(e)]
valid_errs = [e for e in errors if e is not None and not np.isnan(e)]

axes[0].loglog(valid_h, valid_errs, 'o-', color='steelblue', linewidth=2)
axes[0].set_xlabel('h (step size)')
axes[0].set_ylabel('|estimate - true limit|')
axes[0].set_title('Error in numerical limit\nBreakdown at very small h')
axes[0].invert_xaxis()

# Optimal h zone illustration
h_range = np.logspace(-15, 0, 300)
# truncation error decreases as h->0, roundoff error increases
truncation = h_range          # O(h) for first-order approx
roundoff   = 1e-16 / h_range  # machine epsilon / h
total = truncation + roundoff

axes[1].loglog(h_range, truncation, '--', color='orange', label='Truncation error O(h)')
axes[1].loglog(h_range, roundoff,   '--', color='red',    label='Roundoff error ε/h')
axes[1].loglog(h_range, total,      '-',  color='steelblue', linewidth=2, label='Total error')
axes[1].axvline(1e-8, color='green', linestyle=':', label='Optimal h ≈ √ε')
axes[1].set_xlabel('h')
axes[1].set_ylabel('Error')
axes[1].set_title('Truncation vs roundoff tradeoff\n(optimal h ≈ √machine_epsilon ≈ 1e-8)')
axes[1].legend(fontsize=9)

plt.tight_layout()
plt.show()

def check_limit(f, a, steps=15, tol=1e-6):
    """Check whether lim_{x->a} f(x) exists by comparing left and right limits."""
    h_vals = 1.0 / (2.0 ** np.arange(steps))
    
    right_vals = np.array([f(a + h) for h in h_vals])
    left_vals  = np.array([f(a - h) for h in h_vals])
    
    # Use median of last 5 values to reduce noise
    right_limit = np.median(right_vals[-5:])
    left_limit  = np.median(left_vals[-5:])
    
    exists = abs(right_limit - left_limit) < tol
    return {
        'right_limit': right_limit,
        'left_limit':  left_limit,
        'limit_exists': exists,
        'limit_value':  (right_limit + left_limit) / 2 if exists else None
    }


cases = [
    ('sin(x)/x at x=0',       lambda x: np.sin(x)/x, 0.0),
    ('|x|/x at x=0 (sign)',   lambda x: np.abs(x)/x, 0.0),
    ('x^2 at x=2',            lambda x: x**2,         2.0),
    ('1/x at x=0',            lambda x: 1/x,           0.0),
]

for name, f, a in cases:
    try:
        result = check_limit(f, a)
        print(f"{name}")
        print(f"  Right limit: {result['right_limit']:.6f}")
        print(f"  Left limit:  {result['left_limit']:.6f}")
        print(f"  Limit exists: {result['limit_exists']}" +
              (f", value: {result['limit_value']:.6f}" if result['limit_exists'] else ''))
        print()
    except Exception as e:
        print(f"{name}: Error — {e}\n")

# Softmax temperature limiting behavior
def softmax(x, temperature=1.0):
    z = x / temperature
    z = z - np.max(z)  # numerical stability
    e = np.exp(z)
    return e / e.sum()

logits = np.array([2.0, 1.0, 0.5, -0.5])
temps  = [10.0, 1.0, 0.5, 0.1, 0.01]

print('Softmax output as temperature T → 0')
print(f'Logits: {logits}')
print()
print(f'{"T":>8}  {"p[0]":>8}  {"p[1]":>8}  {"p[2]":>8}  {"p[3]":>8}')
print('-' * 50)
for T in temps:
    p = softmax(logits, T)
    print(f'{T:>8.2f}  {p[0]:>8.4f}  {p[1]:>8.4f}  {p[2]:>8.4f}  {p[3]:>8.4f}')

print()
print('As T→0, the distribution concentrates on the maximum logit (index 0).')
print('In the limit T→0: softmax → one-hot at argmax.')
print('As T→∞: softmax → uniform distribution.')