第 2 章 極限的嚴格化：柯西、魏爾斯特拉斯與 ε-δ#

import sympy as sp
from fractions import Fraction

# ---- 狄德金切割：把「實數」定義成一個有理數的下界集合 ----
def sqrt2_cut(p):
    """sqrt(2) 的切割：精確有理下界 q，滿足 q <= sqrt(2) < q + 10^-p"""
    scale = 10**p
    q = sp.floor(sp.sqrt(2 * scale**2))
    return Fraction(int(q), scale)

print("狄德金切割：sqrt(2) 的下界集合")
for p in range(1, 6):
    q = sqrt2_cut(p)
    print(f"    p = {p}   下界 q = {str(q):>22}   q^2 <= 2 ? {q*q <= 2}   誤差 = {float(sp.sqrt(2)) - float(q):.3e}")
print()

# 切割的運算：加法是集合的 Minkowski 和
A = [Fraction(1, 1), Fraction(3, 2)]
B = [Fraction(1, 2), Fraction(2, 3)]
print("切割的加法：A + B =", [str(v) for v in sorted({a + b for a in A for b in B})])
print()

# ---- e 的有理數逼近：柯西數列在切割的觀點下收斂到「一個實數」 ----
print("e 的有理數逼近（切割的上半部）:")
for N in [3, 6, 12]:
    approx = sum((Fraction(1, int(sp.factorial(k))) for k in range(0, N + 1)), Fraction(0))
    print(f"    N = {N:>3}  部分和 = {float(approx):.15f}   真值 e = {float(sp.E):.15f}   誤差 = {abs(float(approx)-float(sp.E)):.3e}")
print()

# 柯西準則：相鄰項的差距必須趨於 0（切割定義的可用性）
seq_e = [(1 + 1 / sp.Integer(m))**m for m in range(1, 60)]
print("柯西數列 (1+1/n)^n 的相鄰項差距（必須趨於 0）：")
for m in [10, 20, 30, 40]:
    print(f"    n = {m:>3}   a_(n+1) - a_n = {sp.N(seq_e[m] - seq_e[m-1], 8)}")

狄德金切割：sqrt(2) 的下界集合
    p = 1   下界 q =                    7/5   q^2 <= 2 ? True   誤差 = 1.421e-02
    p = 2   下界 q =                141/100   q^2 <= 2 ? True   誤差 = 4.214e-03
    p = 3   下界 q =                707/500   q^2 <= 2 ? True   誤差 = 2.136e-04
    p = 4   下界 q =              7071/5000   q^2 <= 2 ? True   誤差 = 1.356e-05
    p = 5   下界 q =          141421/100000   q^2 <= 2 ? True   誤差 = 3.562e-06

切割的加法：A + B = ['3/2', '5/3', '2', '13/6']

e 的有理數逼近（切割的上半部）:
    N =   3  部分和 = 2.666666666666667   真值 e = 2.718281828459045   誤差 = 5.162e-02
    N =   6  部分和 = 2.718055555555555   真值 e = 2.718281828459045   誤差 = 2.263e-04
    N =  12  部分和 = 2.718281828286169   真值 e = 2.718281828459045   誤差 = 1.729e-10

柯西數列 (1+1/n)^n 的相鄰項差距（必須趨於 0）：
    n =  10   a_(n+1) - a_n = 0.010456552
    n =  20   a_(n+1) - a_n = 0.0029655088
    n =  30   a_(n+1) - a_n = 0.0013775300
    n =  40   a_(n+1) - a_n = 0.00079250915

import numpy as np
import sympy as sp

# ---- Heine 的有限覆蓋引理：緊緻集的有限子覆蓋 ----
# 故事：Heine 1824 年用「有限覆蓋」把閉區間上連續函數的最大值定理
# 從歐拉時代的「直覺」變成可檢查的有限步驟。

def greedy_subcover(cover, a, b):
    """從開覆蓋中貪心地抽出有限子覆蓋（覆蓋緊緻集時總能成功）"""
    chosen, pts = [], np.linspace(a, b, 20001)
    uncovered = list(pts)
    changed = True
    while uncovered and changed:
        changed = False
        for (lo, hi) in cover:
            hit = [p for p in uncovered if lo <= p <= hi]
            if hit:
                chosen.append((lo, hi))
                uncovered = [p for p in uncovered if not (lo <= p <= hi)]
                changed = True
                if not uncovered:
                    break
    return chosen, uncovered

# 一個開覆蓋：半徑隨位置變化的區間
cover = sorted((k * 0.05 - (0.06 + 0.02 * abs(k * 0.05)),
                k * 0.05 + (0.06 + 0.02 * abs(k * 0.05)))
               for k in range(-2, 60))

fin, left = greedy_subcover(cover, 0.0, 1.0)
print("Heine 有限覆蓋引理的數值版本")
print("  開覆蓋中區間的總數 :", len(cover))
print("  貪心抽出的有限子覆蓋 :", len(fin), "個")
print("  仍未覆蓋的點數 :", len(left))
gaps = [fin[i][1] - fin[i + 1][0] for i in range(len(fin) - 1)]
print("  相鄰子覆蓋區間的最大『缺口』 =", max(gaps))
print()

# ---- Bolzano–Weierstrass：任何有界數列都有收斂子列 ----
# 定理的標準證明是「對半分割 + 鴿籠」：把 [lo,hi] 對半切，
# 永遠保留含項較多的那一半，區間長度 -> 0，其中無窮多項收斂到一點。
def bw_nested(a, rounds=30):
    lo, hi = float(a.min()), float(a.max())
    keep = a
    history = []
    for _ in range(rounds):
        mid = (lo + hi) / 2
        left, right = keep[keep <= mid], keep[keep > mid]
        if len(left) >= len(right):
            keep, hi = left, mid
        else:
            keep, lo = right, mid
        history.append((lo, hi, len(keep)))
    return keep, (lo + hi) / 2, history

# 一個有聚點的數列：偶數項恰為 0，奇數項趨近 0
a = np.concatenate([np.zeros(150), 1.0 / np.arange(1, 150)])
rng = np.random.default_rng(20261106)
a = np.concatenate([a, rng.uniform(-1, 1, 100)])

keep, L, hist = bw_nested(a)
print("Bolzano–Weierstrass 定理（對半分割法）")
print("  原數列長度 :", len(a))
print("  第 0 步區間長度 :", hist[0][1] - hist[0][0], " 留存項數 :", hist[0][2])
for r in [1, 5, 11]:
    lo, hi, cnt = hist[r]
    print(f"    第 {r+1:>2} 步: 區間 = [{lo:.3e}, {hi:.3e}], 長度 = {hi-lo:.3e}, 留存項數 = {cnt}")
print("  極限 L =", L)
print("  留存子列的極差 =", np.ptp(keep))
print("  留存子列是否收斂到 L ?", np.allclose(keep, L, atol=1e-6))
print()

# ---- 單調有界定理：單調 + 有界 => 收斂 ----
n = sp.symbols('n', integer=True, positive=True)
b_n = 1 - 1 / n
print("單調有界定理的符號示範：a_n = 1 - 1/n")
print("     a_(n+1) - a_n =", sp.simplify(b_n.subs(n, n + 1) - b_n))
print("     1 - a_n =", sp.simplify(1 - b_n))
print("     limit =", sp.limit_seq(b_n, n))

Heine 有限覆蓋引理的數值版本
  開覆蓋中區間的總數 : 62
  貪心抽出的有限子覆蓋 : 21 個
  仍未覆蓋的點數 : 0
  相鄰子覆蓋區間的最大『缺口』 = 0.10699999999999987

Bolzano–Weierstrass 定理（對半分割法）
  原數列長度 : 399
  第 0 步區間長度 : 0.9821254977259336  留存項數 : 284
    第  2 步: 區間 = [-4.732e-01, 1.787e-02], 長度 = 4.911e-01, 留存項數 = 263
    第  6 步: 區間 = [-1.282e-02, 1.787e-02], 長度 = 3.069e-02, 留存項數 = 244
    第 12 步: 區間 = [-3.485e-04, 1.310e-04], 長度 = 4.796e-04, 留存項數 = 150
  極限 L = -4.987592959796879e-10
  留存子列的極差 = 0.0
  留存子列是否收斂到 L ? True

單調有界定理的符號示範：a_n = 1 - 1/n
     a_(n+1) - a_n = 1/(n*(n + 1))
     1 - a_n = 1/n
     limit = 1

import numpy as np
import sympy as sp

# ---- 魏爾斯特拉斯 1860 年代的核心：連續函數的極限要『一致收斂』才保連續 ----
x, n = sp.symbols('x n', real=True)

# (a) 逐點收斂：連續函數列的極限函數可以不連續
print("逐點收斂的陷阱：f_n(x) = x^n 在 [0,1] 上")
for xx in [sp.Rational(1, 2), sp.Rational(3, 4), sp.Integer(1)]:
    val = sp.limit(xx**n, n, sp.oo)
    print(f"    逐點極限在 x = {xx} :", val)
print("    極限函數 f(x) = 0 (x<1), f(1) = 1，在 x=1 不連續")
print()

# (b) 一致收斂：極限函數必連續
print("一致收斂的正確版本：g_n(x) = (1 - x/n)^n 在 [0,1] 上")
for xx in [sp.Rational(1, 2), sp.Integer(1)]:
    val = sp.limit((1 - xx / n)**n, n, sp.oo)
    print(f"    逐點極限在 x = {xx} :", val)
print("    極限函數 = exp(-x)，連續")
print("    檢驗 exp(-x) 在 x -> 1^+ :", sp.limit(sp.exp(-x), x, 1, dir='+'))
print()

# ---- 一致收斂的『速率』：用數值檢查 sup |g_n - e^{-x}| -> 0 ----
xs = np.linspace(0, 1, 2001)
def sup_err(N):
    gn = (1 - xs / N) ** N
    return float(np.max(np.abs(gn - np.exp(-xs))))
print("  一致收斂的定量檢查：sup_{x in [0,1]} |g_n(x) - e^{-x}|")
for N in [1, 2, 5, 10, 50]:
    print(f"    n = {N:>4}  sup 誤差 = {sup_err(N):.6e}")
print("  -> 隨 n 趨於 0，這就是『一致』的數值意義")
print()

# ---- 連續函數在緊緻集上必一致連續（Heine–Cantor）----
eps = 1e-6
def delta_needed(x0, eps):
    # |x^2 - x0^2| = |x-x0|(x+x0) <= |x-x0|(2*x0+delta)
    # 取 delta 使 delta*(2*x0+delta) < eps
    return (-2 * x0 + np.sqrt(4 * x0 ** 2 + 4 * eps)) / 2
xs2 = np.linspace(0, 2, 1001)
deltas = np.array([delta_needed(xx, eps) for xx in xs2])
print("Heine–Cantor：連續函數在緊緻集上必一致連續")
print("    f(x) = x^2, eps =", eps)
print("    各點所需的 delta 的最小值 =", deltas.min())
print("    各點所需的 delta 的最大值 =", deltas.max())
print("    比值 max/min =", deltas.max() / deltas.min())
print("    -> delta 可以『統一』選取，這就是一致連續")

逐點收斂的陷阱：f_n(x) = x^n 在 [0,1] 上
    逐點極限在 x = 1/2 : 0
    逐點極限在 x = 3/4 : 0
    逐點極限在 x = 1 : 1
    極限函數 f(x) = 0 (x<1), f(1) = 1，在 x=1 不連續

一致收斂的正確版本：g_n(x) = (1 - x/n)^n 在 [0,1] 上
    逐點極限在 x = 1/2 : exp(-1/2)
    逐點極限在 x = 1 : exp(-1)
    極限函數 = exp(-x)，連續
    檢驗 exp(-x) 在 x -> 1^+ : exp(-1)

  一致收斂的定量檢查：sup_{x in [0,1]} |g_n(x) - e^{-x}|
    n =    1  sup 誤差 = 3.678794e-01
    n =    2  sup 誤差 = 1.178794e-01
    n =    5  sup 誤差 = 4.019944e-02
    n =   10  sup 誤差 = 1.920100e-02
    n =   50  sup 誤差 = 3.709761e-03
  -> 隨 n 趨於 0，這就是『一致』的數值意義

Heine–Cantor：連續函數在緊緻集上必一致連續
    f(x) = x^2, eps = 1e-06
    各點所需的 delta 的最小值 = 2.4999998426977754e-07
    各點所需的 delta 的最大值 = 0.001
    比值 max/min = 4000.0002516835752
    -> delta 可以『統一』選取，這就是一致連續

import numpy as np
import sympy as sp

# ---- 魏爾斯特拉斯的『不等式的演算術』：delta 是被『構造』出來的，不是被『找到』的 ----
x = sp.symbols('x', real=True)

def eps_delta_check(f, a, eps, delta, trials=8000, seed=20261106):
    """對候選 delta 做隨機抽查：確認 |f(x)-f(a)| < eps 在 0<|x-a|<delta 時成立"""
    rng = np.random.default_rng(seed)
    xs = a + rng.uniform(-delta, delta, trials)
    xs = xs[xs != a]
    vals = np.array([float(sp.N(f.subs(x, xv))) for xv in xs])
    base = float(sp.N(f.subs(x, a)))
    worst = float(np.max(np.abs(vals - base)))
    return worst, worst < eps

print("魏爾斯特拉斯的『不等式的演算術』：delta 是被『構造』出來的，不是被『找到』的")
print()
cases = [
    ("f(x) = x^2,    a = 2", x**2,      sp.Integer(2), lambda e: e / 5),
    ("f(x) = sin(x), a = 0", sp.sin(x), sp.Integer(0), lambda e: e),
    ("f(x) = 1/x,    a = 1", 1 / x,      sp.Integer(1), lambda e: e / 2),
    ("f(x) = exp(x), a = 0", sp.exp(x), sp.Integer(0), lambda e: sp.log(1 + e)),
]
for name, f, a, rule in cases:
    print(name)
    for eps in [sp.Rational(1, 10), sp.Rational(1, 10**6)]:
        delta = rule(eps)
        worst, ok = eps_delta_check(f, a, eps, delta)
        print(f"    eps = {str(eps):>12}   delta = {float(delta):.3e}   抽查最大誤差 = {worst:.3e}   成立 ? {ok}")
    print()

# ---- 對稱化技巧：處理 x = 1/x 這種『不能直接取 delta』的情形 ----
print("對稱化技巧：|1/x - 1/a| = |x-a| / (|x| |a|)")
print("  a = 1，若 |x-1| < 1/2 則 |x| > 1/2，故 |1/x - 1| < 2|x-1|，取 delta = eps/2")
eps = sp.Rational(1, 1000)
delta = eps / 2
worst, ok = eps_delta_check(1 / x, sp.Integer(1), eps, delta)
print("    eps =", eps, "  delta =", delta, "  抽查最大誤差 =", worst, "  成立 ?", ok)
print()

# ---- SymPy 自動化的不等式求解器可以幫我們『找出』delta ----
print("SymPy 的不等式求解器：直接解出『什麼 x 會讓 |f(x)-f(a)| < eps』")
y = sp.symbols('y', real=True)
print("  解 (y^2 - 4) < eps :", sp.solve_univariate_inequality((y**2 - 4) < eps, y, relational=False))
print("  解 (y**2 - 4) > eps :", sp.solve_univariate_inequality((y**2 - 4) > eps, y, relational=False))
print("  -> 這是『算不等式』，正是柯西 -> 魏爾斯特拉斯這條線的核心技術")

魏爾斯特拉斯的『不等式的演算術』：delta 是被『構造』出來的，不是被『找到』的

f(x) = x^2,    a = 2
    eps =         1/10   delta = 2.000e-02   抽查最大誤差 = 8.040e-02   成立 ? True
    eps =    1/1000000   delta = 2.000e-07   抽查最大誤差 = 8.000e-07   成立 ? True

f(x) = sin(x), a = 0
    eps =         1/10   delta = 1.000e-01   抽查最大誤差 = 9.983e-02   成立 ? True
    eps =    1/1000000   delta = 1.000e-06   抽查最大誤差 = 1.000e-06   成立 ? True

f(x) = 1/x,    a = 1
    eps =         1/10   delta = 5.000e-02   抽查最大誤差 = 5.261e-02   成立 ? True
    eps =    1/1000000   delta = 5.000e-07   抽查最大誤差 = 5.000e-07   成立 ? True

f(x) = exp(x), a = 0
    eps =         1/10   delta = 9.531e-02   抽查最大誤差 = 1.000e-01   成立 ? True
    eps =    1/1000000   delta = 1.000e-06   抽查最大誤差 = 1.000e-06   成立 ? True

對稱化技巧：|1/x - 1/a| = |x-a| / (|x| |a|)
  a = 1，若 |x-1| < 1/2 則 |x| > 1/2，故 |1/x - 1| < 2|x-1|，取 delta = eps/2
    eps = 1/1000   delta = 1/2000   抽查最大誤差 = 0.0005000502104153259   成立 ? True

SymPy 的不等式求解器：直接解出『什麼 x 會讓 |f(x)-f(a)| < eps』
  解 (y^2 - 4) < eps : Interval.open(-sqrt(40010)/100, sqrt(40010)/100)
  解 (y**2 - 4) > eps : Union(Interval.open(-oo, -sqrt(40010)/100), Interval.open(sqrt(40010)/100, oo))
  -> 這是『算不等式』，正是柯西 -> 魏爾斯特拉斯這條線的核心技術