13.2 數值分析的誤差理論與計算機時代的微積分#

import sys
import numpy as np
import sympy as sp

print("=== 1. 浮點的基本事實 ===")
print("float64 的 machine epsilon =", np.finfo(np.float64).eps)
print("float64 的尾數位數 mant_dig =", sys.float_info.mant_dig)
print("符號有理數 sp.Rational(1,3) =", sp.Rational(1, 3))
print("轉成 float 後 =", repr(float(sp.Rational(1, 3))))
print("SymPy 的 evalf 可以任意精度：")
for d in [15, 30, 50, 80]:
    print(f"   sp.N(sp.sqrt(2), {d}) =", sp.N(sp.sqrt(2), d))

print()
print("=== 2. 相對誤差如何被放大 ===")
x = 1.0 + 1e-12
print("x = 1 + 1e-12 存成 float =", repr(x))
print("x - 1 =", repr(x - 1.0), " （相對誤差被放大了 1e12 倍）")
print("相對誤差 =", (x - 1.0) / 1.0)

print()
print("=== 3. 相消誤差（catastrophic cancellation）===")
a, b = 1.0, 1.0 - 1e-17
print("a =", repr(a), " b =", repr(b), " a - b =", repr(a - b))
print("精確的 a - b 應該是 1e-17，但 float 給出", repr(a - b))
print("SymPy 精確版本 a - b =", sp.Rational(1) - (sp.Rational(1) - sp.Rational(1, 10**17)))

print()
print("=== 4. SymPy 用 Rational + evalf 做高精度 ===")
q = sp.Rational(1, 3)
print("Rational(1,3) 的 30 位小數 =", sp.N(q, 30))
print("pi 的 50 位 =", sp.N(sp.pi, 50))
print("e 的 30 位 =", sp.N(sp.E, 30))

=== 1. 浮點的基本事實 ===
float64 的 machine epsilon = 2.220446049250313e-16
float64 的尾數位數 mant_dig = 53
符號有理數 sp.Rational(1,3) = 1/3
轉成 float 後 = 0.3333333333333333
SymPy 的 evalf 可以任意精度：
   sp.N(sp.sqrt(2), 15) = 1.41421356237310
   sp.N(sp.sqrt(2), 30) = 1.41421356237309504880168872421
   sp.N(sp.sqrt(2), 50) = 1.4142135623730950488016887242096980785696718753769
   sp.N(sp.sqrt(2), 80) = 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070

=== 2. 相對誤差如何被放大 ===
x = 1 + 1e-12 存成 float = 1.000000000001
x - 1 = 1.000088900582341e-12  （相對誤差被放大了 1e12 倍）
相對誤差 = 1.000088900582341e-12

=== 3. 相消誤差（catastrophic cancellation）===
a = 1.0  b = 1.0  a - b = 0.0
精確的 a - b 應該是 1e-17，但 float 給出 0.0
SymPy 精確版本 a - b = 1/100000000000000000

=== 4. SymPy 用 Rational + evalf 做高精度 ===
Rational(1,3) 的 30 位小數 = 0.333333333333333333333333333333
pi 的 50 位 = 3.1415926535897932384626433832795028841971693993751
e 的 30 位 = 2.71828182845904523536028747135

import numpy as np
import sympy as sp

# ============================================================
# Wilkinson（1969）：病態線性方程組與條件數
# ============================================================
A = np.array([[1.0, 1.0], [1.0, 1.0 + 1e-10]])
b = np.array([1.0, 1.0 + 1e-10])
x0 = np.linalg.solve(A, b)
print("A =", A.tolist())
print("b =", b.tolist())
print("解（浮點） =", x0)
print("A 的 2-範數條件數 =", np.linalg.cond(A))

print()
print("擾動 b 的第二個分量，看解的相對變化：")
for delta in [1e-16, 1e-14, 1e-12]:
    bp = b.copy()
    bp[1] += delta
    xp = np.linalg.solve(A, bp)
    rel = abs(xp[1] - x0[1]) / abs(x0[1])
    print(f"   擾動 {delta:.0e} -> x2 = {xp[1]:.12f}   相對變化 = {rel:.3e}")

print()
print("用 SymPy 精確求解對照：")
As = sp.Matrix([[1, 1], [1, sp.Rational(10000000001, 10**10)]])
bs = sp.Matrix([1, sp.Rational(10000000001, 10**10)])
xs = As.LUsolve(bs)
print("   精確解 =", [sp.N(v, 20) for v in xs])

A = [[1.0, 1.0], [1.0, 1.0000000001]]
b = [1.0, 1.0000000001]
解（浮點） = [0. 1.]
A 的 2-範數條件數 = 39999965677.40396

擾動 b 的第二個分量，看解的相對變化：
   擾動 1e-16 -> x2 = 1.000000000000   相對變化 = 0.000e+00
   擾動 1e-14 -> x2 = 1.000099920064   相對變化 = 9.992e-05
   擾動 1e-12 -> x2 = 1.010000888178   相對變化 = 1.000e-02

用 SymPy 精確求解對照：
   精確解 = [0, 1.0000000000000000000]

import numpy as np
import sympy as sp

# ============================================================
# (a) 牛頓法 vs 二分法：收斂率
# ============================================================
z = sp.Symbol('z')
root_hi = sp.nsolve(z**3 - 2*z - 5, 2.0, n=60)
root = float(root_hi)
print("高精度參考根（SymPy nsolve, 60 位）=", root_hi)


def f(x): return x**3 - 2*x - 5
def df(x): return 3*x**2 - 2


def newton(f, df, x0, maxit=6):
    x = x0
    hist = []
    for k in range(maxit):
        hist.append((k, x, abs(f(x))))
        x = x - f(x)/df(x)
    return hist


hist = newton(f, df, 2.0)
print()
print(f"{'k':>3} {'x_k':>22} {'|f(x_k)|':>12} {'誤差 e_k':>12} {'e_k/e_{k-1}':>14}")
prev = None
for k, x, r in hist:
    e = abs(x - root)
    ratio = (e/prev) if prev else float('nan')
    print(f"{k:>3} {x:>22.16f} {r:>12.3e} {e:>12.3e} {ratio:>14.3e}")
    prev = e


def bisect(f, a, b, tol=1e-12, maxit=200):
    fa = f(a)
    for k in range(maxit):
        m = 0.5*(a+b)
        fm = f(m)
        if abs(fm) < tol:
            return m, k
        if fa*fm < 0:
            b = m
        else:
            a, fa = m, fm
    return 0.5*(a+b), maxit


xb, kb = bisect(f, 2.0, 3.0)
print()
print(f"二分法：{kb} 次迭代 -> x = {xb!r}")
print(f"   誤差 = {abs(xb-root):.3e}；區間寬 2^-{kb} = {2.0**(-kb):.3e}")

# ============================================================
# (b) Taylor 餘項：截斷誤差（第 5 章回扣）
# ============================================================
x = sp.Symbol('x')
print()
print("用符號 Taylor 展開算 sin(1)：")
for n in [3, 5, 7, 11, 15, 21]:
    approx = sp.series(sp.sin(x), x, 0, n).removeO().subs(x, 1)
    print(f"   n={n:>3} 近似 = {sp.N(approx, 18)}   誤差 = {float(abs(approx - sp.sin(1))):.3e}")

# ============================================================
# (c) 截斷誤差 vs 捨入誤差（第 5.3 節的經典實驗）
# ============================================================
print()
print("用 float 直接累加 sin(1) 的 Taylor 級數：")
s = 0.0
term = 1.0
for k in range(0, 30):
    if k > 0:
        term = term * (-1.0)/((2*k)*(2*k+1))
    s += term
    if k in (2, 4, 6, 8, 12, 16, 20, 25):
        print(f"   取到 k={k:>2} 項   累積和 = {s!r}   誤差 = {abs(s - float(sp.sin(1))):.3e}")

# ============================================================
# (d) 災難性相消：float 版 exp(-30) 的 Taylor 級數
# ============================================================
print()
print("用 float 算 exp(-30) 的 Taylor 級數（項的絕對值遠大於結果）：")
x0 = -30.0
term = 1.0
s = 1.0
for k in range(1, 61):
    term = term * x0 / k
    s += term
    if k in (10, 20, 30, 40, 50, 60):
        print(f"   取到 k={k:>3} 項   累積和 = {s!r}   誤差 = {abs(s - float(sp.exp(-30.0))):.3e}")
print("   真值 =", repr(float(sp.exp(-30))))

高精度參考根（SymPy nsolve, 60 位）= 2.09455148154233

  k                    x_k     |f(x_k)|       誤差 e_k    e_k/e_{k-1}
  0     2.0000000000000000    1.000e+00    9.455e-02            nan
  1     2.1000000000000001    6.100e-02    5.449e-03      5.762e-02
  2     2.0945681211041851    1.857e-04    1.664e-05      3.054e-03
  3     2.0945514816981992    1.740e-09    1.559e-10      9.368e-06
  4     2.0945514815423265    8.882e-16    0.000e+00      0.000e+00
  5     2.0945514815423265    8.882e-16    0.000e+00            nan

二分法：41 次迭代 -> x = 2.0945514815423394
   誤差 = 1.288e-14；區間寬 2^-41 = 4.547e-13

用符號 Taylor 展開算 sin(1)：
   n=  3 近似 = 1.00000000000000000   誤差 = 1.585e-01
   n=  5 近似 = 0.833333333333333333   誤差 = 8.138e-03
   n=  7 近似 = 0.841666666666666667   誤差 = 1.957e-04
   n= 11 近似 = 0.841471009700176367   誤差 = 2.489e-08
   n= 15 近似 = 0.841470984808658420   誤差 = 7.619e-13
   n= 21 近似 = 0.841470984807896507   誤差 = 1.953e-20

用 float 直接累加 sin(1) 的 Taylor 級數：
   取到 k= 2 項   累積和 = 0.8416666666666667   誤差 = 1.957e-04
   取到 k= 4 項   累積和 = 0.8414710097001764   誤差 = 2.489e-08
   取到 k= 6 項   累積和 = 0.8414709848086585   誤差 = 7.619e-13
   取到 k= 8 項   累積和 = 0.8414709848078965   誤差 = 0.000e+00
   取到 k=12 項   累積和 = 0.8414709848078965   誤差 = 0.000e+00
   取到 k=16 項   累積和 = 0.8414709848078965   誤差 = 0.000e+00
   取到 k=20 項   累積和 = 0.8414709848078965   誤差 = 0.000e+00
   取到 k=25 項   累積和 = 0.8414709848078965   誤差 = 0.000e+00

用 float 算 exp(-30) 的 Taylor 級數（項的絕對值遠大於結果）：
   取到 k= 10 項   累積和 = 121254849.57142858   誤差 = 1.213e+08
   取到 k= 20 項   累積和 = 85291712210.03218   誤差 = 8.529e+10
   取到 k= 30 項   累積和 = 384842612535.9078   誤差 = 3.848e+11
   取到 k= 40 項   累積和 = 63336540117.26651   誤差 = 6.334e+10
   取到 k= 50 項   累積和 = 878229229.2779531   誤差 = 8.782e+08
   取到 k= 60 項   累積和 = 1685584.3003382375   誤差 = 1.686e+06
   真值 = 9.357622968840175e-14

import numpy as np
import sympy as sp
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

x = sp.Symbol('x')
f_sym = sp.exp(x)
A, B = 0.0, 1.0
exact = sp.integrate(f_sym, (x, A, B))
exact_f = float(exact)
print("區間 [0,1] 上 exp(x) 的精確值 =", exact, " =", sp.N(exact, 20))


def f(t): return np.exp(t)


def trapezoid(f, a, b, n):
    xs = np.linspace(a, b, n+1)
    ys = f(xs)
    return (b-a)/n * (0.5*ys[0] + ys[1:-1].sum() + 0.5*ys[-1])


def romberg(f, a, b, kmax=6):
    """R[k][0] = 2^k 段的梯形法；R[k][j] = k+1 階 Richardson 外推。"""
    R = np.zeros((kmax, kmax))
    R[0, 0] = trapezoid(f, a, b, 1)
    for k in range(1, kmax):
        h = (b-a)/2**k
        R[k, 0] = 0.5*R[k-1, 0] + h*f(a + h*np.arange(1, 2**k, 2)).sum()
        for j in range(1, k+1):
            R[k, j] = R[k, j-1] + (R[k, j-1] - R[k-1, j-1])/(4**j - 1)
    return R


R = romberg(f, A, B, 6)
print()
print("Romberg 表：")
for i in range(6):
    print("   ", "  ".join(f"{R[i, j]:.15f}" for j in range(i+1)))
print()
print("對角線 R[k][k] 的誤差：")
for k in range(6):
    print(f"   R[{k}][{k}] = {R[k, k]:.16f}   誤差 = {abs(R[k, k]-exact_f):.3e}")


def gauss_legendre(f, a, b, n):
    nodes, weights = np.polynomial.legendre.leggauss(n)
    half = (b-a)/2
    mid = (a+b)/2
    return half*np.sum(weights*f(mid + half*nodes))


print()
print("Gauss-Legendre 求積：")
gl_err = []
for n in [1, 2, 4, 8, 12, 16]:
    v = gauss_legendre(f, A, B, n)
    gl_err.append(abs(v-exact_f))
    print(f"   n={n:>3}  值 = {v:.16f}  誤差 = {abs(v-exact_f):.3e}")

print()
print("複合梯形法（對照）：")
tp_err = []
for n in [2, 4, 8, 16, 32, 64]:
    v = trapezoid(f, A, B, n)
    tp_err.append(abs(v-exact_f))
    print(f"   n={n:>4}  值 = {v:.16f}  誤差 = {abs(v-exact_f):.3e}")

# 畫圖：誤差隨求值次數的變化
fig, ax = plt.subplots(figsize=(7, 4.5))
ns = np.array([2, 4, 8, 16, 32, 64])
ax.loglog(ns, tp_err, "o-", label="composite trapezoid  (O(h^2))")
ax.loglog([2, 4, 8, 16], [gl_err[0], gl_err[1], gl_err[3], gl_err[5]],
           "s-", label="Gauss-Legendre  (n points -> O(h^{2n-1}))")
rb = [abs(R[k, k]-exact_f) for k in range(6)]
ax.loglog([2**k for k in range(6)], rb, "^-", label="Romberg diagonal")
ax.set_xlabel("number of function evaluations / panels")
ax.set_ylabel("absolute error")
ax.set_title("Quadrature error on int_0^1 exp(x) dx")
ax.legend()
ax.grid(alpha=0.3, which="both")
plt.tight_layout()
plt.savefig("/tmp/ch13_quad.png", dpi=120)
print("saved /tmp/ch13_quad.png")

區間 [0,1] 上 exp(x) 的精確值 = 1.71828182845905  = 1.7182818284590450908

Romberg 表：
    1.859140914229523
    1.753931092464825  1.718861151876593
    1.727221904557517  1.718318841921747  1.718282687924757
    1.720518592164302  1.718284154699897  1.718281842218440  1.718281828794530
    1.718841128579994  1.718281974051892  1.718281828675358  1.718281828460389  1.718281828459078
    1.718421660316327  1.718281837561772  1.718281828462430  1.718281828459051  1.718281828459045  1.718281828459045

對角線 R[k][k] 的誤差：
   R[0][0] = 1.8591409142295225   誤差 = 1.409e-01
   R[1][1] = 1.7188611518765928   誤差 = 5.793e-04
   R[2][2] = 1.7182826879247572   誤差 = 8.595e-07
   R[3][3] = 1.7182818287945303   誤差 = 3.355e-10
   R[4][4] = 1.7182818284590784   誤差 = 3.331e-14
   R[5][5] = 1.7182818284590453   誤差 = 2.220e-16

Gauss-Legendre 求積：
   n=  1  值 = 1.6487212707001282  誤差 = 6.956e-02
   n=  2  值 = 1.7178963780075041  誤差 = 3.855e-04
   n=  4  值 = 1.7182818275260778  誤差 = 9.330e-10
   n=  8  值 = 1.7182818284590451  誤差 = 0.000e+00
   n= 12  值 = 1.7182818284590453  誤差 = 2.220e-16
   n= 16  值 = 1.7182818284590455  誤差 = 4.441e-16

複合梯形法（對照）：
   n=   2  值 = 1.7539310924648253  誤差 = 3.565e-02
   n=   4  值 = 1.7272219045575166  誤差 = 8.940e-03
   n=   8  值 = 1.7205185921643018  誤差 = 2.237e-03
   n=  16  值 = 1.7188411285799945  誤差 = 5.593e-04
   n=  32  值 = 1.7184216603163267  誤差 = 1.398e-04
   n=  64  值 = 1.7183167868500937  誤差 = 3.496e-05
saved /tmp/ch13_quad.png

import math
import numpy as np
import sympy as sp

# ============================================================
# Kahan 補償求和（1965）與 Neumaier 改進
# ============================================================
def naive_sum(vals):
    s = 0.0
    for v in vals:
        s += v
    return s


def kahan_sum(vals):
    """Kahan（1965）：把每一步的捨入誤差存起來，下一步減掉。"""
    s = 0.0
    c = 0.0
    for v in vals:
        y = v - c
        t = s + y
        c = (t - s) - y
        s = t
    return s


def neumaier_sum(vals):
    """Neumaier（Kahan–Babuška–Neumaier）：對大小不一的加數更穩健。"""
    s = 0.0
    c = 0.0
    for v in vals:
        t = s + v
        if abs(s) >= abs(v):
            c += (s - t) + v
        else:
            c += (v - t) + s
        s = t
    return s + c


# --- 案例 1：教科書式的相消 ---
vals = [1.0, 1e16, 1.0, -1e16]
print("vals =", vals)
print("  naive   =", naive_sum(vals))
print("  kahan   =", kahan_sum(vals))
print("  neumaier=", neumaier_sum(vals))
print("  SymPy 精確 = 2")

# --- 案例 2：0.1 加十次 ---
vals2 = [0.1]*10
print()
print("0.1 加 10 次：")
print("  naive   =", repr(naive_sum(vals2)))
print("  kahan   =", repr(kahan_sum(vals2)))
print("  neumaier=", repr(neumaier_sum(vals2)))
print("  SymPy 精確 =", sp.Rational(1, 10)*10)

# --- 案例 3：等比級數，累積誤差隨 N 成長 ---
print()
N = 200000
r = 1.0000001
g = [1e-8 * r**k for k in range(N)]
print("N =", N, " 的等比級數 sum 1e-8 * r^k")
print("  naive   =", naive_sum(g))
print("  kahan   =", kahan_sum(g))
print("  neumaier=", neumaier_sum(g))
print("  math.fsum =", math.fsum(g))
print("  閉式精確值 =", 1e-8*(1 - r**N)/(1 - r))
print()
print("累積誤差 |naive - neumaier| 隨 N 成長：")
for n in [10**3, 10**4, 10**5, N]:
    vv = g[:n]
    print(f"   N={n:>8}  |naive - neumaier| = {abs(naive_sum(vv)-neumaier_sum(vv)):.3e}")

vals = [1.0, 1e+16, 1.0, -1e+16]
  naive   = 0.0
  kahan   = 0.0
  neumaier= 2.0
  SymPy 精確 = 2

0.1 加 10 次：
  naive   = 0.9999999999999999
  kahan   = 1.0
  neumaier= 1.0
  SymPy 精確 = 1

N = 200000  的等比級數 sum 1e-8 * r^k
  naive   = 0.002020133900667315
  kahan   = 0.002020133900667288
  neumaier= 0.002020133900667288
  math.fsum = 0.002020133900667288
  閉式精確值 = 0.002020133900667298

累積誤差 |naive - neumaier| 隨 N 成長：
   N=    1000  |naive - neumaier| = 1.016e-20
   N=   10000  |naive - neumaier| = 8.132e-20
   N=  100000  |naive - neumaier| = 2.602e-18
   N=  200000  |naive - neumaier| = 2.689e-17

完全精確                                              完全不精確
  |<---------------------- SymPy ---------------------->|
   Rational, Integer, Symbol, Add, Mul, Pow ...     evalf 到 n 位
        (精確算術，零誤差)                        (高精度數值近似)

import numpy as np
import sympy as sp

# ============================================================
# 符號 vs 數值：分工的實際示範
# ============================================================
n = sp.Symbol('n', integer=True, positive=True)

# --- 1) 符號端：精確封閉形式 ---
s1 = sp.summation(1/n**2, (n, 1, sp.oo))
print("sum 1/n^2  =", s1, "=", sp.N(s1, 20))
s2 = sp.summation(sp.Rational(1, 2)**n, (n, 1, sp.oo))
print("sum (1/2)^n =", s2)

# --- 2) 符號端：精確有理數運算 ---
print()
print("Rational(1,3) + Rational(1,6) =", sp.Rational(1, 3) + sp.Rational(1, 6), "（精確，不是浮點）")
print("   =", sp.N(sp.Rational(1, 3) + sp.Rational(1, 6), 20), "（轉成小數）")

# --- 3) 數值端：直接浮點累加會怎樣 ---
print()
vals = [1.0/n**2 for n in range(1, 200001)]
s_float = 0.0
for v in vals:
    s_float += v
s_exact = float(s1)
print("sum_{n=1}^{200000} 1/n^2 的浮點累加 =", repr(s_float))
print("sum_{n=1}^{inf} 1/n^2（符號精確）  =", s_exact)
print("差 =", abs(s_float - s_exact))

# --- 4) 符號端：積分也能精確到封閉形式 ---
x = sp.Symbol('x')
I = sp.integrate(sp.exp(-x**2), (x, -sp.oo, sp.oo))
print()
print("int exp(-x^2) dx（整個實軸）= ", I, "=", sp.N(I, 20))
print("這在數值上需要無窮多點；符號運算一步到位。")

# --- 5) 混合：符號求導 + 高精度數值代入 ---
print()
f = x**2 * sp.exp(-x)
df = sp.diff(f, x)
print("f'(x) =", df)
for d in [15, 30, 50]:
    print(f"   f'(1.7) 到 {d} 位 =", sp.N(df.subs(x, sp.Rational(17, 10)), d))

sum 1/n^2  = pi**2/6 = 1.6449340668482264365
sum (1/2)^n = 1

Rational(1,3) + Rational(1,6) = 1/2 （精確，不是浮點）
   = 0.50000000000000000000 （轉成小數）

sum_{n=1}^{200000} 1/n^2 的浮點累加 = 1.6449290668607415
sum_{n=1}^{inf} 1/n^2（符號精確）  = 1.6449340668482264
差 = 4.999987484932689e-06

int exp(-x^2) dx（整個實軸）=  sqrt(pi) = 1.7724538509055160273
這在數值上需要無窮多點；符號運算一步到位。

f'(x) = -x**2*exp(-x) + 2*x*exp(-x)
   f'(1.7) 到 15 位 = 0.0931685972668947
   f'(1.7) 到 30 位 = 0.0931685972668946716141894272571
   f'(1.7) 到 50 位 = 0.093168597266894671614189427257059341242538795175326