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

import sympy as sp

# ---- 夾擠定理（Squeeze / Sandwich）----
x, n = sp.symbols('x n', real=True, positive=True)

print("=== 夾擠定理 ===")
print("目標：lim(x->0) sin(x)/x = 1")
print("  夾擠：cos(x) <= sin(x)/x <= 1  (0 < x < pi/2)")
print("  cos(x) 的極限 =", sp.limit(sp.cos(x), x, 0))
print("  1 的極限      =", sp.limit(sp.Integer(1), x, 0))
print("  所以夾擠定理給出 lim sin(x)/x = 1")
print()

t = sp.symbols('t', positive=True)
print("  SymPy 驗證不等式 sin(t) - t*cos(t) >= 0 的級數展開：")
print("    t - sin(t)        =", sp.series(t - sp.sin(t), t, 0, 5))
print("    sin(t) - t*cos(t) =", sp.series(sp.sin(t) - t * sp.cos(t), t, 0, 5))
print()

print("  數值夾擠驗證：")
for xv in [sp.Rational(1, 1000), sp.Rational(1, 100), sp.Rational(1, 10)]:
    lo = sp.N(sp.cos(xv), 12)
    mid = sp.N(sp.sin(xv) / xv, 12)
    hi = sp.Integer(1)
    print(f"    x = {str(xv):>6}: {lo} <= {mid} <= 1  ? {lo <= mid <= 1}")
print()

# ---- 單調有界定理 ----
print("=== 單調有界定理 ===")
m = sp.symbols('m', integer=True, positive=True)
a_m = 1 - 1 / m
print("  a_m = 1 - 1/m")
print("    a_{m+1} - a_m =", sp.simplify(a_m.subs(m, m + 1) - a_m), " (> 0 所以遞增)")
print("    1 - a_m =", sp.simplify(1 - a_m), " (> 0 所以 a_m < 1 有界)")
print("    limit_seq =", sp.limit_seq(a_m, m))

=== 夾擠定理 ===
目標：lim(x->0) sin(x)/x = 1
  夾擠：cos(x) <= sin(x)/x <= 1  (0 < x < pi/2)
  cos(x) 的極限 = 1
  1 的極限      = 1
  所以夾擠定理給出 lim sin(x)/x = 1

  SymPy 驗證不等式 sin(t) - t*cos(t) >= 0 的級數展開：
    t - sin(t)        = t**3/6 + O(t**5)
    sin(t) - t*cos(t) = t**3/3 + O(t**5)

  數值夾擠驗證：
    x = 1/1000: 0.999999500000 <= 0.999999833333 <= 1  ? True
    x =  1/100: 0.999950000417 <= 0.999983333417 <= 1  ? True
    x =   1/10: 0.995004165278 <= 0.998334166468 <= 1  ? True

=== 單調有界定理 ===
  a_m = 1 - 1/m
    a_{m+1} - a_m = 1/(m*(m + 1))  (> 0 所以遞增)
    1 - a_m = 1/m  (> 0 所以 a_m < 1 有界)
    limit_seq = 1

import sympy as sp

m = sp.symbols('m', integer=True, positive=True)

# ---- Stolz-Cesaro 的等價形式 ----
# 定理：若 (b_n) 嚴格遞增且 b_n -> +inf，且 lim (a_{n+1}-a_n)/(b_{n+1}-b_n) = L 存在，
#       則 lim a_n/b_n = L。
# 應用 1：(1+1/n)^n
#   取 a_k = ln(1+1/k), b_k = 1/k（皆遞減，等價於對 -a_k, -b_k 使用遞增型）。
A = sp.log(1 + 1 / m)
B = 1 / m
dA = sp.simplify(A.subs(m, m + 1) - A)
dB = sp.simplify(B.subs(m, m + 1) - B)
print("Stolz-Cesaro 應用 1：證 lim n*ln(1+1/n) = 1")
print("  a_{n+1} - a_n =", dA)
print("  b_{n+1} - b_n =", dB)
ratio = sp.simplify(dA / dB)
print("  差商 =", ratio)
print("  lim 差商 =", sp.limit(ratio, m, sp.oo))
print("  -> lim n*ln(1+1/n) = 1, 故 lim (1+1/n)^n = exp(1) =", sp.exp(1))
print()

# 應用 2：ln n / n -> 0 (用於證 lim n^(1/n) = 1)
print("Stolz-Cesaro 應用 2：證 lim ln(n)/n = 0")
dA2 = sp.simplify(sp.log(m + 1) - sp.log(m))
dB2 = sp.simplify((m + 1) - m)
print("  a_{n+1} - a_n =", dA2, " b_{n+1} - b_n =", dB2)
print("  lim 差商 =", sp.limit(dA2 / dB2, m, sp.oo))
print("  -> lim ln(n)/n = 0, 故 lim n^(1/n) = exp(0) = 1")
print()

# 應用 3：2^n / n! -> 0
print("Stolz-Cesaro 應用 3：證 lim 2^n/n! = 0")
a3 = 2**m
b3 = sp.factorial(m)
dA3 = sp.simplify(a3.subs(m, m + 1) - a3)
dB3 = sp.simplify(b3.subs(m, m + 1) - b3)
ratio3 = sp.simplify(dA3 / dB3)
print("  a_{n+1} - a_n =", dA3)
print("  b_{n+1} - b_n =", dB3)
print("  差商 =", ratio3)
print("  lim 差商 =", sp.limit(ratio3, m, sp.oo))
print("  -> lim 2^n/n! = 0")
print()

# 應用 4：H_n / ln n -> 1（Stolz 的招牌例子）
print("Stolz-Cesaro 應用 4：證 lim H_n/ln(n) = 1")
H = sp.harmonic(m)
dA4 = sp.simplify(H.subs(m, m + 1) - H)
dB4 = sp.simplify(sp.log(m + 1) - sp.log(m))
print("  a_{n+1} - a_n =", dA4)
print("  b_{n+1} - b_n =", dB4)
print("  lim 差商 =", sp.limit(dA4 / dB4, m, sp.oo))

Stolz-Cesaro 應用 1：證 lim n*ln(1+1/n) = 1
  a_{n+1} - a_n = log(m*(m + 2)/(m + 1)**2)
  b_{n+1} - b_n = -1/(m*(m + 1))
  差商 = -log((m*(m + 2)/(m + 1)**2)**(m*(m + 1)))
  lim 差商 = 1
  -> lim n*ln(1+1/n) = 1, 故 lim (1+1/n)^n = exp(1) = E

Stolz-Cesaro 應用 2：證 lim ln(n)/n = 0
  a_{n+1} - a_n = log((m + 1)/m)  b_{n+1} - b_n = 1
  lim 差商 = 0
  -> lim ln(n)/n = 0, 故 lim n^(1/n) = exp(0) = 1

Stolz-Cesaro 應用 3：證 lim 2^n/n! = 0
  a_{n+1} - a_n = 2**m
  b_{n+1} - b_n = m*factorial(m)
  差商 = 2**m/(m*factorial(m))
  lim 差商 = 0
  -> lim 2^n/n! = 0

Stolz-Cesaro 應用 4：證 lim H_n/ln(n) = 1
  a_{n+1} - a_n = -harmonic(m) + harmonic(m + 1)
  b_{n+1} - b_n = log((m + 1)/m)
  lim 差商 = 1

import sympy as sp

n = sp.symbols('n', integer=True, positive=True)
x = sp.symbols('x', real=True)

print("=== 四個招牌極限的 SymPy 驗證 ===")
targets = [
    ("lim n*sin(1/n)", n * sp.sin(1 / n), 1),
    ("lim (1+1/n)^n",  (1 + 1 / n)**n,     sp.E),
    ("lim 2^n/n!",     2**n / sp.factorial(n), 0),
    ("lim n^(1/n)",    n**(1 / n),        1),
]
for name, expr, expect in targets:
    got = sp.limit(expr, n, sp.oo)
    print(f"  {name:<22} = {str(got):<8}  期望 {str(expect):<6}  相符 ? {sp.simplify(got - expect) == 0}")
print()

print("=== 對數化簡 ===")
for name, expr in [("(1+1/n)^n", (1 + 1 / n)**n),
                   ("n^(1/n)",    n**(1 / n)),
                   ("2^n/n!",     2**n / sp.factorial(n))]:
    lg = sp.limit(sp.log(expr), n, sp.oo)
    print(f"  {name:<14}:  ln 後的極限 = {str(lg):<6} -> exp 回來 = {sp.exp(lg)}")
print()

print("=== dir 引數：從哪一邊逼近 ===")
print("  lim(x->0) 1/x        dir='+'  :", sp.limit(1 / x, x, 0, dir='+'))
print("  lim(x->0) 1/x        dir='-'  :", sp.limit(1 / x, x, 0, dir='-'))
print("  lim(x->0) 1/x        dir='+-':", sp.limit(1 / x, x, 0, dir='+-'))
print("  lim(x->0) |x|/x     dir='+'  :", sp.limit(sp.Abs(x) / x, x, 0, dir='+'))
print("  lim(x->0) |x|/x     dir='-'  :", sp.limit(sp.Abs(x) / x, x, 0, dir='-'))
try:
    sp.limit(sp.Abs(x) / x, x, 0, dir='+-')
except ValueError as e:
    print("  lim(x->0) |x|/x     dir='+-':  ValueError ->", e)
print()

print("=== sp.Limit：未求值的極限物件 ===")
L = sp.Limit(sp.sin(x) / x, x, 0)
print("  L =", L)
print("  L.doit() =", L.doit())
print("  L 的 args :", L.args)
print()

print("=== function_range：函數在區間上的值域 ===")
from sympy.calculus.util import function_range
for f, iv in [(x**2, sp.Interval(-3, 3)),
              (x**2 - 2 * x, sp.Interval(0, 4))]:
    print(f"  f(x) = {str(f):<16} 在 {str(iv):<18} 的值域 =", function_range(f, x, iv))
print()

print("=== AccumBounds：SymPy 的『誠實失敗』 ===")
A = sp.limit(sp.exp(-x**2) * sp.sin(1 / x), x, 0)
print("  lim(x->0) e^{-x^2} sin(1/x) =", A)
print("  它的上下界 =", A.args)
B = sp.limit(sp.sin(x) / x, x, 0)
print("  能算出來時給精確值 :", B, " / 類別 :", type(B).__name__)
print("  算不出來時給弱結論 :", A, " / 類別 :", type(A).__name__)

=== 四個招牌極限的 SymPy 驗證 ===
  lim n*sin(1/n)         = 1         期望 1       相符 ? True
  lim (1+1/n)^n          = E         期望 E       相符 ? True
  lim 2^n/n!             = 0         期望 0       相符 ? True
  lim n^(1/n)            = 1         期望 1       相符 ? True

=== 對數化簡 ===
  (1+1/n)^n     :  ln 後的極限 = 1      -> exp 回來 = E
  n^(1/n)       :  ln 後的極限 = 0      -> exp 回來 = 1
  2^n/n!        :  ln 後的極限 = -oo    -> exp 回來 = 0

=== dir 引數：從哪一邊逼近 ===
  lim(x->0) 1/x        dir='+'  : oo
  lim(x->0) 1/x        dir='-'  : -oo
  lim(x->0) 1/x        dir='+-': zoo
  lim(x->0) |x|/x     dir='+'  : 1
  lim(x->0) |x|/x     dir='-'  : -1
  lim(x->0) |x|/x     dir='+-':  ValueError -> The limit does not exist since left hand limit = -1 and right hand limit = 1

=== sp.Limit：未求值的極限物件 ===
  L = Limit(sin(x)/x, x, 0, dir='+')
  L.doit() = 1
  L 的 args : (sin(x)/x, x, 0, +)

=== function_range：函數在區間上的值域 ===
  f(x) = x**2             在 Interval(-3, 3)    的值域 = Interval(0, 9)
  f(x) = x**2 - 2*x       在 Interval(0, 4)     的值域 = Interval(-1, 8)

=== AccumBounds：SymPy 的『誠實失敗』 ===
  lim(x->0) e^{-x^2} sin(1/x) = AccumBounds(-1, 1)
  它的上下界 = (-1, 1)
  能算出來時給精確值 : 1  / 類別 : One
  算不出來時給弱結論 : AccumBounds(-1, 1)  / 類別 : AccumulationBounds

The limit does not exist since left hand limit = -1 and right hand limit = 1

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

# ---- 用數值與隨機取樣尋找反例的工作流 ----
# 定理候選：「lim(x->0) sin(x)/x = 1」
# 我們不去證明它，而是先努力把它『推翻』：隨機取樣，看有沒有反例。
rng = np.random.default_rng(115)
x = sp.symbols('x', real=True)

# (1) 隨機取樣：在 (0, 1) 上撒 2000 個點，看 |f(x)-1| 會不會超過某個門檻
xs = rng.uniform(1e-6, 1.0, 2000)
errs = np.array([abs(np.sin(xv) / xv - 1.0) for xv in xs])
order = np.argsort(xs)
half = len(xs) // 2
print("(1) 隨機取樣尋找反例：lim(x->0) sin(x)/x = 1")
print("    取樣點數 :", len(xs))
print("    最大誤差 :", errs.max(), "  對應的 x =", xs[np.argmax(errs)])
print("    最小 x   :", xs.min(), "（最靠近 0 的取樣點）")
print("    最小一半 x 的最大誤差 :", errs[order[:half]].max())
print("    較大一半 x 的最大誤差 :", errs[order[half:]].max())
print("    -> 靠近 0 的一批誤差遠小於遠離 0 的一批，方向正確")
print()

# (2) 更聰明的取樣：對 x 用對數刻度，才能逼近 0
logxs = np.logspace(-12, 0, 6)
print("(2) 對數刻度取樣（逼近 0 的正確方式）")
for xv in logxs:
    print(f"    x = {xv:>10.3e}   sin(x)/x - 1 = {np.sin(xv)/xv - 1:+.6e}   x^2/6 = {xv**2/6:.6e}")
print("    -> 誤差與 x^2/6 同階，這預告了誤差是 O(x^2)")
print()

# (3) 用 SymPy 反過來驗證這個預測：級數展開
print("(3) SymPy 級數展開確認誤差量級")
print("    sin(x)/x - 1 =", sp.series(sp.sin(x) / x - 1, x, 0, 7))
print()

# (4) 估計極限值，並與猜的係數比對：候選命題是錯的
print("(4) 用取樣估計極限值，並與猜的係數比對")
print("    候選命題：lim(x->0) (e^x - 1 - x)/x^2 = 1/6 （猜錯了）")
xs3 = rng.uniform(1e-4, 1e-3, 500)
vals = np.array([(np.exp(xv) - 1 - xv) / xv**2 for xv in xs3])
print("    隨機取樣估計出的極限值 ≈", round(float(vals.mean()), 6), " （而不是 1/6 =", 1/6, "）")
print("    -> 差異很大，說明候選命題是假的")
print("    SymPy 級數揭穿真相：")
print("      ", sp.series((sp.exp(x) - 1 - x) / x**2, x, 0, 5))
print("    SymPy 直接算出的極限 =", sp.limit((sp.exp(x) - 1 - x) / x**2, x, 0))
print()

# (5) 找到一個明確的反例：把『對所有 x>0 成立』的強命題推翻
print("(5) 推翻強命題：『|sin(x)| <= x^2 對所有 x > 0 成立』")
viol_x = None
for _ in range(2000):
    xv = rng.uniform(1e-6, 2.0)
    if abs(np.sin(xv)) > xv**2:
        viol_x = xv
        break
print("    隨機取樣找到的反例 x =", viol_x)
print("    精確驗證：取 x = 1/2")
print("      |sin(1/2)| =", sp.Abs(sp.sin(sp.Rational(1, 2))))
print("      (1/2)^2    =", sp.Rational(1, 4))
print("      反例成立 ?", sp.Abs(sp.sin(sp.Rational(1, 2))) > sp.Rational(1, 4))
print("    -> 反例的真正用處不是推翻結論，而是提醒我們：")
print("       定理總要帶上『適用範圍』，而找出那個範圍正是嚴格化的意義")
print()

# (6) 把取樣結果畫出來：雙對數座標確認誤差量級
fig, ax = plt.subplots(figsize=(7, 4.5))
xs4 = np.logspace(-10, -1, 400)
ax.loglog(xs4, np.abs(np.sin(xs4) / xs4 - 1), "b-", label="|sin(x)/x - 1|")
ax.loglog(xs4, xs4**2 / 6, "k--", label="x^2 / 6")
ax.loglog(xs4, xs4**2 / 7, "k:", label="x^2 / 7")
ax.set_xlabel("x")
ax.set_ylabel("error")
ax.set_title("Error of sin(x)/x near 0 (log-log scale)")
ax.legend()
plt.tight_layout()
plt.savefig("/tmp/sin_over_x_error.png", dpi=120)
print("    已存 /tmp/sin_over_x_error.png")
print("    兩條參考線夾住實測曲線 -> 誤差確實是 Theta(x^2)")

(1) 隨機取樣尋找反例：lim(x->0) sin(x)/x = 1
    取樣點數 : 2000
    最大誤差 : 0.1578013418560268   對應的 x = 0.9975815110891472
    最小 x   : 1.5983754111375778e-05 （最靠近 0 的取樣點）
    最小一半 x 的最大誤差 : 0.04178966839166842
    較大一半 x 的最大誤差 : 0.1578013418560268
    -> 靠近 0 的一批誤差遠小於遠離 0 的一批，方向正確

(2) 對數刻度取樣（逼近 0 的正確方式）
    x =  1.000e-12   sin(x)/x - 1 = +0.000000e+00   x^2/6 = 1.666667e-25
    x =  2.512e-10   sin(x)/x - 1 = +0.000000e+00   x^2/6 = 1.051596e-20
    x =  6.310e-08   sin(x)/x - 1 = -6.661338e-16   x^2/6 = 6.635120e-16
    x =  1.585e-05   sin(x)/x - 1 = -4.186473e-11   x^2/6 = 4.186477e-11
    x =  3.981e-03   sin(x)/x - 1 = -2.641487e-06   x^2/6 = 2.641489e-06
    x =  1.000e+00   sin(x)/x - 1 = -1.585290e-01   x^2/6 = 1.666667e-01
    -> 誤差與 x^2/6 同階，這預告了誤差是 O(x^2)

(3) SymPy 級數展開確認誤差量級
    sin(x)/x - 1 = -x**2/6 + x**4/120 - x**6/5040 + O(x**7)

(4) 用取樣估計極限值，並與猜的係數比對
    候選命題：lim(x->0) (e^x - 1 - x)/x^2 = 1/6 （猜錯了）
    隨機取樣估計出的極限值 ≈ 0.50009  （而不是 1/6 = 0.16666666666666666 ）
    -> 差異很大，說明候選命題是假的
    SymPy 級數揭穿真相：
       1/2 + x/6 + x**2/24 + x**3/120 + x**4/720 + O(x**5)
    SymPy 直接算出的極限 = 1/2

(5) 推翻強命題：『|sin(x)| <= x^2 對所有 x > 0 成立』
    隨機取樣找到的反例 x = 0.8393914180193854
    精確驗證：取 x = 1/2
      |sin(1/2)| = sin(1/2)
      (1/2)^2    = 1/4
      反例成立 ? True
    -> 反例的真正用處不是推翻結論，而是提醒我們：
       定理總要帶上『適用範圍』，而找出那個範圍正是嚴格化的意義

    已存 /tmp/sin_over_x_error.png
    兩條參考線夾住實測曲線 -> 誤差確實是 Theta(x^2)