12.1 勒貝格測度與積分#

import sympy as sp
from fractions import Fraction
import math

# ---- Jordan 外含量：對有限個有理數點，用任意小的總長覆蓋 ----
eps = Fraction(1, 100)
pts = [Fraction(k, 37) for k in range(1, 8)]
w = Fraction(eps, 2 * len(pts))
cover = [(p - w / 2, p + w / 2) for p in pts]
total = sum((r - l) for l, r in cover)
print('覆蓋的點數 =', len(pts))
print('每個區間寬度 =', w)
print('總長度 =', total, '  eps =', eps, '  總長 < eps :', total < eps)

# ---- Lebesgue 外測度：用可數多個等長區間覆蓋 [0,1]，取所有覆蓋的下確界 ----
print()
for m in [10, 100, 1000]:
    length = Fraction(1, m)
    n_needed = math.ceil(1 / length)
    print(f'用長度 1/{m} 的區間覆蓋 [0,1]，需要 {n_needed} 個，總長 = {n_needed*length}')

# ---- Jordan 含量只能有限可加，無法可數可加 ----
print()
print('Jordan 含量(Q ∩ [0,1])   = 0')
print('Jordan 含量([0,1] \\ Q)   = 0')
print('但 Jordan 含量([0,1])     = 1')
print('=> [0,1] 分割成兩個含量皆為 0 的集合，Jordan 含量對「可數分割」無法保持可加性')

覆蓋的點數 = 7
每個區間寬度 = 1/1400
總長度 = 1/200   eps = 1/100   總長 < eps : True

用長度 1/10 的區間覆蓋 [0,1]，需要 10 個，總長 = 1
用長度 1/100 的區間覆蓋 [0,1]，需要 100 個，總長 = 1
用長度 1/1000 的區間覆蓋 [0,1]，需要 1000 個，總長 = 1

Jordan 含量(Q ∩ [0,1])   = 0
Jordan 含量([0,1] \ Q)   = 0
但 Jordan 含量([0,1])     = 1
=> [0,1] 分割成兩個含量皆為 0 的集合，Jordan 含量對「可數分割」無法保持可加性

import sympy as sp

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

# Vitali 集合 V：V 與 V+q (q 有理) 互不相交，且 [0,1] 被 {V+q} 覆蓋
print('情形一：若 m(V) = m > 0，則可數多個平移的測度總和 =', sp.summation(m, (n, 1, sp.oo)))
print('        但 [0,1] 的測度只有 1，而平移的並集覆蓋 [0,1]，矛盾')
print('情形二：若 m(V) = 0，則平移的測度總和 = 0，但並集應等於整個 [0,1]，矛盾')
print('=> V 不可測：兩種情形都不可能')

# ---- 良序化的算術核心：Vitali 集合的存在正是良序化（= 選擇公理）的直接後果 ----
print()
print('可數多個 1 的和 =', sp.summation(sp.Integer(1), (n, 1, sp.oo)))
print('可數多個 0 的和 =', sp.summation(sp.Integer(0), (n, 1, sp.oo)))

情形一：若 m(V) = m > 0，則可數多個平移的測度總和 = oo
        但 [0,1] 的測度只有 1，而平移的並集覆蓋 [0,1]，矛盾
情形二：若 m(V) = 0，則平移的測度總和 = 0，但並集應等於整個 [0,1]，矛盾
=> V 不可測：兩種情形都不可能

可數多個 1 的和 = oo
可數多個 0 的和 = 0

import sympy as sp

# ---- 可數集的外測度 = 0：逐點用越來越小的區間覆蓋 ----
eps = sp.Rational(1, 100)
pts = [sp.Rational(k, 97) for k in range(1, 6)]
widths = [eps / 2**(k + 1) for k in range(len(pts))]
print('被覆蓋的點   =', pts)
print('各覆蓋區間長度 =', widths)
print('總長 =', sum(widths), '  <  eps =', eps, ' :', sum(widths) < eps)
print('=> 可數集的外測度 m* = 0（因為 eps 可任意小）')

# ---- Caratheodory 條件對 E = [0,1/2] 成立 ----
print()
a, b = sp.Rational(0), sp.Rational(1)
e = sp.Rational(1, 2)
print('m*(A∩E) =', e - a, '   m*(A\\E) =', b - e,
      '   總和 =', (e - a) + (b - e), '   m*(A) =', b - a)
print('Caratheodory 條件成立：', (e - a) + (b - e) == b - a)

# ---- Lebesgue 測度在可測集上滿足可數可加 ----
print()
print('m([0,1/2]) + m([1/2,1]) =', sp.Rational(1, 2) + sp.Rational(1, 2),
      ' = m([0,1]) =', sp.Integer(1))
print('m([0,1] ∩ Q) = 0（Q 可數）  但  m([0,1]) = 1  ->  「零測度集」不是「空集」')

# ---- 連續函數：黎曼和的極限 = 勒貝格積分 = sp.integrate 的結果 ----
print()
k = sp.symbols('k', integer=True, positive=True)
Nn = sp.symbols('Nn', integer=True, positive=True)
x = sp.symbols('x', real=True)
R = sp.summation((k / Nn)**2 / Nn, (k, 1, Nn))
print('Riemann 和 =', R)
print('極限       =', sp.limit(R, Nn, sp.oo))
print('勒貝格積分 =', sp.integrate(x**2, (x, 0, 1)))

被覆蓋的點   = [1/97, 2/97, 3/97, 4/97, 5/97]
各覆蓋區間長度 = [1/200, 1/400, 1/800, 1/1600, 1/3200]
總長 = 31/3200   <  eps = 1/100  : True
=> 可數集的外測度 m* = 0（因為 eps 可任意小）

m*(A∩E) = 1/2    m*(A\E) = 1/2    總和 = 1    m*(A) = 1
Caratheodory 條件成立： True

m([0,1/2]) + m([1/2,1]) = 1  = m([0,1]) = 1
m([0,1] ∩ Q) = 0（Q 可數）  但  m([0,1]) = 1  ->  「零測度集」不是「空集」

Riemann 和 = (Nn**3/3 + Nn**2/2 + Nn/6)/Nn**3
極限       = 1/3
勒貝格積分 = 1/3

import sympy as sp

def dirichlet(xv):
    """xv 為 sympy 實數：有理回傳 1，無理回傳 0"""
    return sp.Integer(1) if xv.is_rational else sp.Integer(0)

def riemann_sum(f, a, b, N, sampler):
    h = (b - a) / N
    return sp.simplify(sum(f(sampler(a + k * h)) for k in range(N)) * h)

for N in [4, 8, 16]:
    r = riemann_sum(dirichlet, sp.Rational(0), sp.Rational(1), N, lambda t: t)
    i = riemann_sum(dirichlet, sp.Rational(0), sp.Rational(1), N, lambda t: t + sp.sqrt(2) / 7)
    print(f'N={N:3d}  有理採樣的和 = {r}   無理採樣的和 = {i}')

print('黎曼積分：上和 = 1，下和 = 0，永遠不相等 -> 黎曼不可積')
print('勒貝格積分：1*m(Q∩[0,1]) + 0*m([0,1]\\Q) = 1*0 + 0*1 =', sp.Integer(0))

N=  4  有理採樣的和 = 1   無理採樣的和 = 0
N=  8  有理採樣的和 = 1   無理採樣的和 = 0
N= 16  有理採樣的和 = 1   無理採樣的和 = 0
黎曼積分：上和 = 1，下和 = 0，永遠不相等 -> 黎曼不可積
勒貝格積分：1*m(Q∩[0,1]) + 0*m([0,1]\Q) = 1*0 + 0*1 = 0

from fractions import Fraction

def svc_intervals(n):
    """Smith-Volterra-Cantor 集：第 n 步之後剩下的閉區間"""
    ivs = [(Fraction(0), Fraction(1))]
    for k in range(1, n + 1):
        gap = Fraction(1, 4 ** k)
        new = []
        for (a, b) in ivs:
            mid = (a + b) / 2
            new.append((a, mid - gap / 2))
            new.append((mid + gap / 2, b))
        ivs = new
    return ivs

def darboux(ivs, N):
    """在均勻 N 分割上計算下和與上和（Darboux 逼近）"""
    lower = 0.0
    upper = 0.0
    h = 1.0 / N
    for i in range(N):
        lo, hi = Fraction(i, N), Fraction(i + 1, N)
        hit = False
        inside = True
        for (a, b) in ivs:
            if a < hi and lo < b:
                hit = True
                if not (a <= lo and hi <= b):
                    inside = False
        if hit:
            upper += h
        if hit and inside:
            lower += h
    return lower, upper

print('Smith-Volterra-Cantor 集的剩餘總長（= Lebesgue 測度的近似）:')
for n in [1, 3, 8]:
    tot = sum(b - a for a, b in svc_intervals(n))
    print(f'  n={n}: 總長 = {tot} ({float(tot):.6f})   -> 極限 1/2')

print()
for N in [200, 800]:
    lo, up = darboux(svc_intervals(10), N)
    print(f'均勻 {N} 分割：下和 = {lo:.4f}   上和 = {up:.4f}   振幅 = {up - lo:.4f}')
print('振幅不趨於 0 -> 1_S 不是黎曼可積的；但 Lebesgue 積分 = m(S) = 1/2')

Smith-Volterra-Cantor 集的剩餘總長（= Lebesgue 測度的近似）:
  n=1: 總長 = 3/4 (0.750000)   -> 極限 1/2
  n=3: 總長 = 9/16 (0.562500)   -> 極限 1/2
  n=8: 總長 = 257/512 (0.501953)   -> 極限 1/2

均勻 200 分割：下和 = 0.0000   上和 = 0.6000   振幅 = 0.6000
均勻 800 分割：下和 = 0.0000   上和 = 0.5500   振幅 = 0.5500
振幅不趨於 0 -> 1_S 不是黎曼可積的；但 Lebesgue 積分 = m(S) = 1/2

import sympy as sp

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

# (1) Fatou：移動尖峰 f_n = n * 1_{(0,1/n)}，積分恆為 1，但極限函數幾乎處處為 0
I_move = sp.integrate(n, (x, 0, 1 / n))
print('int f_n =', I_move, ' -> liminf =', sp.limit(I_move, n, sp.oo))
print('int liminf f_n (a.e. 為 0) =', sp.Integer(0), ' <= 1 : Fatou 的不等式嚴格成立')

# (2) 控制收斂：f_n = min(n, 1/sqrt(x))，被可積的 g = 1/sqrt(x) 控制
I1 = sp.integrate(n, (x, 0, 1 / n**2))
I2 = sp.integrate(x**sp.Rational(-1, 2), (x, 1 / n**2, 1))
tot = sp.simplify(I1 + I2)
print()
print('int f_n =', tot, '  -> lim =', sp.limit(tot, n, sp.oo))
print('int g   =', sp.integrate(x**sp.Rational(-1, 2), (x, 0, 1)), '（g 可積）')

# (3) 控制收斂：f_n = x^n (1-x)，被 x(1-x) 控制
J = sp.simplify(sp.integrate(x**n * (1 - x), (x, 0, 1)))
print()
print('int x^n(1-x) =', J, '  -> lim =', sp.limit(J, n, sp.oo))
print('lim f_n = 0 (a.e.)，int lim = 0')

# (4) 單調收斂（遞減版本）：f_n = 1 - x^n
K = sp.simplify(sp.integrate(1 - x**n, (x, 0, 1)))
print()
print('int (1 - x^n) =', K, '  -> lim =', sp.limit(K, n, sp.oo))
print('lim f_n = 1 (a.e.)，int lim = 1  -> 等號成立')

int f_n = 1  -> liminf = 1
int liminf f_n (a.e. 為 0) = 0  <= 1 : Fatou 的不等式嚴格成立

int f_n = 2 - 1/n   -> lim = 2
int g   = 2 （g 可積）

int x^n(1-x) = 1/(n**2 + 3*n + 2)   -> lim = 0
lim f_n = 0 (a.e.)，int lim = 0

int (1 - x^n) = n/(n + 1)   -> lim = 1
lim f_n = 1 (a.e.)，int lim = 1  -> 等號成立

import sympy as sp
import mpmath as mp

x = sp.symbols('x', real=True)
y = sp.symbols('y', real=True)

# (1) 一般的 Fubini：雙重積分先對 x 再對 y
print('int_0^1 int_0^1 (x + y) dx dy =',
      sp.integrate(sp.integrate(x + y, (x, 0, 1)), (y, 0, 1)))

# (2) Tonelli：非負函數，兩個順序必然相等
f = sp.exp(x * y)
ix = sp.integrate(sp.integrate(f, (x, 0, 1)), (y, 0, 1))
iy = sp.integrate(sp.integrate(f, (y, 0, 1)), (x, 0, 1))
print('int int e^(xy) 先 x 後 y =', sp.simplify(ix))
print('int int e^(xy) 先 y 後 x =', sp.simplify(iy))
print('兩者相等 :', sp.simplify(ix - iy) == 0)

# (3) 對角線 {y = x} 的截面測度恆為 0，所以它的二維測度為 0
print()
print('對角線 {y = x} 的二維測度 =', sp.Integer(0), '   而單位正方形的測度 =', sp.Integer(1))

# (4) 沒有絕對可積性時，兩個迭代積分會打架
g = (x**2 - y**2) / (x**2 + y**2)**2
gx = sp.simplify(sp.integrate(sp.integrate(g, (x, 0, 1)), (y, 0, 1)))
gy = sp.simplify(sp.integrate(sp.integrate(g, (y, 0, 1)), (x, 0, 1)))
print()
print('先對 x 積分再對 y :', gx)
print('先對 y 積分再對 x :', gy)
print('兩者相差 =', sp.simplify(gy - gx))

# (5) |g| 的絕對可積性失敗：截去原點附近的正方形後，總量呈對數發散
mp.mp.dps = 20
def num_abs(delta, eps):
    h = lambda a, b: abs((a**2 - b**2) / (a**2 + b**2)**2)
    return mp.quad(h, [delta, eps], [delta, eps])

print()
for delta in [1e-2, 1e-4, 1e-6, 1e-8]:
    val = num_abs(delta, 0.1)
    print(f'截去 [0,{delta}] 後的 int int |g| =', mp.nstr(val, 10),
          '   (pi/4)*log(1/delta) =', mp.nstr(mp.pi / 4 * mp.log(1 / delta), 10))
print('=> |g| 不可積（對數發散），Fubini 的假設不成立，兩個順序自然打架')

int_0^1 int_0^1 (x + y) dx dy = 1
int int e^(xy) 先 x 後 y = -EulerGamma + Ei(1)
int int e^(xy) 先 y 後 x = -EulerGamma + Ei(1)
兩者相等 : True

對角線 {y = x} 的二維測度 = 0    而單位正方形的測度 = 1

先對 x 積分再對 y : -pi/4
先對 y 積分再對 x : pi/4
兩者相差 = pi/2

截去 [0,0.01] 後的 int int |g| = 0.931104909    (pi/4)*log(1/delta) = 3.616892206
截去 [0,0.0001] 後的 int int |g| = 5.33846241    (pi/4)*log(1/delta) = 7.233784412
截去 [0,1e-06] 後的 int int |g| = 9.939943989    (pi/4)*log(1/delta) = 10.85067662
截去 [0,1e-08] 後的 int int |g| = 14.54121853    (pi/4)*log(1/delta) = 14.46756882
=> |g| 不可積（對數發散），Fubini 的假設不成立，兩個順序自然打架