ch220 — Approximation

import numpy as np
import matplotlib.pyplot as plt
from math import factorial

# ── Approximating sin(x) to machine precision ──────────────────────────────
# How many Taylor terms does it take?

def sin_taylor(x, n_terms):
    result = 0.0
    for k in range(n_terms):
        result += ((-1)**k * x**(2*k+1)) / factorial(2*k+1)
    return result

x_test = np.pi / 3  # 60 degrees
exact = np.sin(x_test)

errors = [abs(sin_taylor(x_test, n) - exact) for n in range(1, 12)]
print(f"sin(pi/3) exact: {exact:.15f}")
print()
print(f"{'Terms':>6} | {'Approximation':>18} | {'Error':>12}")
print("-" * 42)
for n in range(1, 12):
    approx = sin_taylor(x_test, n)
    err = abs(approx - exact)
    print(f"{n:>6} | {approx:>18.15f} | {err:>12.2e}")

# Demonstrate truncation vs round-off error tradeoff in numerical diff
# Central difference: (f(x+h) - f(x-h)) / 2h
# Truncation error: O(h^2)  -- smaller h = less truncation
# Round-off error: O(eps/h) -- smaller h = more round-off
# Optimal h ≈ eps^(1/3) for central difference

x0 = 1.2
exact_deriv = np.cos(x0)  # d/dx sin(x) = cos(x)

h_vals = np.logspace(-15, 0, 200)
errors = []
for h in h_vals:
    approx = (np.sin(x0 + h) - np.sin(x0 - h)) / (2*h)
    errors.append(abs(approx - exact_deriv))

plt.figure(figsize=(9, 5))
plt.loglog(h_vals, errors, 'b', lw=1.5)
plt.xlabel('Step size h (log scale)'); plt.ylabel('|error| (log scale)')
plt.title('Truncation vs Round-off: Optimal h for Central Difference')
plt.grid(True, which='both', alpha=0.3)
opt_h = np.finfo(float).eps ** (1/3)
plt.axvline(opt_h, color='red', ls='--', label=f'Optimal h ≈ eps^(1/3) = {opt_h:.1e}')
plt.legend()
plt.tight_layout(); plt.savefig('ch220_tradeoff.png', dpi=100); plt.show()
print(f"Machine epsilon: {np.finfo(float).eps:.2e}")
print(f"Optimal h: {opt_h:.2e}")

# Approximate a noisy signal with polynomials of different degrees
np.random.seed(42)
x_data = np.linspace(0, 2*np.pi, 20)
y_data = np.sin(x_data) + np.random.normal(0, 0.15, 20)

x_plot = np.linspace(0, 2*np.pi, 300)

fig, axes = plt.subplots(1, 3, figsize=(13, 4))
for ax, deg in zip(axes, [2, 5, 15]):
    coeffs = np.polyfit(x_data, y_data, deg)
    y_fit = np.polyval(coeffs, x_plot)
    ax.plot(x_data, y_data, 'ko', ms=5, label='data')
    ax.plot(x_plot, np.sin(x_plot), 'gray', ls='--', label='true sin(x)')
    ax.plot(x_plot, y_fit, 'r', lw=2, label=f'degree {deg}')
    ax.set_ylim(-2.5, 2.5); ax.legend(fontsize=8); ax.grid(True, alpha=0.3)
    ax.set_title(f'Polynomial degree {deg}')
plt.suptitle('Polynomial Approximation: Underfitting, Good Fit, Overfitting')
plt.tight_layout(); plt.savefig('ch220_poly_approx.png', dpi=100); plt.show()