7.3 用 SymPy 做偏導數、重積分與座標變換#

import sympy as sp

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

f = x**2*y - x*sp.sin(y) + sp.exp(x*y)

print('f_x =', sp.diff(f, x))
print('f_y =', sp.diff(f, y))

# 梯度（SymPy 1.x 沒有頂層的 sp.gradient；用 Jacobian 取轉置）
grad = sp.Matrix([f]).jacobian([x, y]).T
print('gradient =')
print(grad)

# Hessian：二階偏導數矩陣
H = sp.hessian(f, (x, y))
print('hessian =')
print(H)

# 混合偏導數相等（Schwarz / Clairaut 定理）
print('f_xy =', sp.simplify(sp.diff(f, x, y)))
print('f_yx =', sp.simplify(sp.diff(f, y, x)))
print('相等？', sp.simplify(sp.diff(f, x, y) - sp.diff(f, y, x)) == 0)

# Hessian 恆為對稱矩陣（Clairaut 定理的矩陣版本）
print('H - H^T =')
print(H - H.T)

# Hessian 的特徵值：決定函數在該點的凹凸與極值型態
q = x**2 + 2*x*y + 3*y**2
Hq = sp.hessian(q, (x, y))
print('二次函數的 Hessian =')
print(Hq)
print('特徵值 =', Hq.eigenvals())
print('det =', Hq.det())
x, y, z = sp.symbols('x y z', real=True)

def grad(f, vars_):
    """回傳 f 對 vars_ 的梯度（列向量 Matrix）。"""
    return sp.Matrix([f]).jacobian(list(vars_)).T

def hess(f, vars_):
    """回傳 f 的 Hessian 矩陣。"""
    return sp.hessian(f, tuple(vars_))

f = x**3*y - sp.sin(x*y)

# 傳入變數元組
print(sp.hessian(f, (x, y)))

# 舊版曾支援的字串寫法 'x,y' 在 SymPy 1.14 已被移除，會拋 ValueError
try:
    sp.hessian(f, 'x,y')
except ValueError as e:
    print('字串寫法已移除:', e)

# 手動取出所有二階偏導數（字典風格，適合程式化處理）
second = {(a, b): sp.diff(f, a, b) for a in (x, y) for b in (x, y)}
print('f_xx =', second[(x, x)])
print('f_xy =', second[(x, y)])
print('f_yx =', second[(y, x)])
print('f_yy =', second[(y, y)])
print('f_xy == f_yx ?', sp.simplify(second[(x, y)] - second[(y, x)]) == 0)

# 三變數版本的輔助函式
g = x**2*sp.sin(y) + y*z**2
print('grad =')
print(grad(g, (x, y, z)))
print('hess =')
print(hess(g, (x, y, z)))
print('grad.f =', sp.simplify((grad(g, (x, y, z)).T * grad(g, (x, y, z)))[0]))

import sympy as sp

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

# 反例：f(x,y) = xy(x^2-y^2)/(x^2+y^2),  f(0,0)=0
f = x*y*(x**2 - y**2)/(x**2 + y**2)

fx = sp.diff(f, x)
fy = sp.diff(f, y)
print('f_x =', fx)
print('f_y =', fy)

# 原點的 f_xy 與 f_yx（用差商定義算，不信任 sp.diff）
fxy0 = sp.limit((fx.subs({y: h, x: 0}) - 0)/h, h, 0)
fyx0 = sp.limit((fy.subs({x: h, y: 0}) - 0)/h, h, 0)
print('f_xy(0,0) =', fxy0)
print('f_yx(0,0) =', fyx0)
print('兩者相等？', sp.simplify(fxy0 - fyx0) == 0)

# 偏導數在原點不連續：沿不同路徑的極限
print('lim_{y->0} f_xy(x,y) (y=0) =', sp.limit(sp.diff(f, x, y).subs({y: 0}), x, 0))
print('lim_{y->0} f_xy(x,y) (x=0) =', sp.limit(sp.diff(f, x, y).subs({x: 0}), y, 0))

# 對照：多項式的混合偏導數恆相等
g = x**3*y - x*y**3 + sp.exp(x*y)
print('g 的 f_xy - f_yx =', sp.simplify(sp.diff(g, x, y) - sp.diff(g, y, x)))

import sympy as sp

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

# 直角座標下的二重積分：矩形區域
f = x**2 + 3*x*y - y**2
val = sp.integrate(f, (y, -2, 2), (x, -1, 1))
print('矩形上的二重積分 =', sp.simplify(val))

# 三角形區域：0 <= x <= 2, 0 <= y <= 2 - x
val2 = sp.integrate(f, (y, 0, 2 - x), (x, 0, 2))
print('三角形上的二重積分 =', sp.simplify(val2))

# 先用 Integral 描述（不積），再 doit()
I = sp.Integral(f, (y, 0, 2 - x), (x, 0, 2))
print('Integral 物件 =', I)
print('doit() =', sp.simplify(I.doit()))

# 圓盤上的二重積分：直角座標（用符號 sqrt 的上下限）
g = x**2 + y**2
disk = sp.integrate(g, (y, -sp.sqrt(1 - x**2), sp.sqrt(1 - x**2)), (x, -1, 1))
print('圓盤上 ∫∫(x^2+y^2) =', sp.simplify(disk))

# --- limits 列表的順序語意（判別式驗證）---
# 區域 0 <= x <= y <= 1：先對 x 積分（0 到 y），再對 y 積分（0 到 1）
h = x + y
v1 = sp.integrate(sp.integrate(h, (x, 0, y)), (y, 0, 1))
print('0<=x<=y<=1 =', sp.simplify(v1))

# 反序：先對 y 積分（x 到 1），再對 x 積分（0 到 1）
v2 = sp.integrate(sp.integrate(h, (y, x, 1)), (x, 0, 1))
print('反序（y 從 x 到 1）= ', sp.simplify(v2))
print('相等？', sp.simplify(v1 - v2) == 0)
print('手算核對 =', sp.Rational(3, 2) * sp.Rational(1, 3))

# 判別式：h2 = x*y，(x,0,y),(y,0,1) 給出 1/8 => x 是內層
h2 = x*y
print('(x,0,y),(y,0,1) =', sp.Integral(h2, (x, 0, y), (y, 0, 1)).doit())
print('(y,0,1),(x,0,y) =', sp.Integral(h2, (y, 0, 1), (x, 0, y)).doit())

# 四分之一圓盤：極座標 vs 直角座標
r, theta = sp.symbols('r theta', real=True)
v3 = sp.integrate(sp.integrate(h.subs({x: r*sp.cos(theta), y: r*sp.sin(theta)})*r, (r, 0, 1)), (theta, 0, sp.pi/2))
print('四分之一圓盤（極座標）= ', sp.simplify(v3))
v4 = sp.integrate(sp.integrate(h, (y, 0, sp.sqrt(1 - x**2))), (x, 0, 1))
print('四分之一圓盤（直角座標）= ', sp.simplify(v4))

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

A = sp.simplify(sp.integrate(sp.integrate(f, (x, 0, y)), (y, 0, 1)))
B = sp.simplify(sp.Integral(f, (x, 0, y), (y, 0, 1)).doit())
print('寫法 A（巢狀）= ', A)
print('寫法 B（Integral）= ', B)
print('A == B ?', sp.simplify(A - B) == 0)

# 真正的陷阱：外層界限裡出現了已被內層積掉的變數
bad = sp.simplify(sp.integrate(sp.integrate(f, (y, 0, 1)), (x, 0, 2-y)))
print('錯誤寫法 =', bad)
print('  -> y 是游離的，這不是任何一個積分的值')

# 正確寫法：先對 x 積分（界限 0 到 2-y，y 還活著），再對 y 積分
good = sp.simplify(sp.integrate(sp.integrate(f, (x, 0, 2-y)), (y, 0, 1)))
print('正確寫法 =', good)

I = sp.Integral(f, (x, 0, 2-y), (y, 0, 1))
print('Integral 結構 =', I)
print('Integral.doit() =', sp.simplify(I.doit()))

import sympy as sp

def change_variables(f, new_x, new_y, old_vars, new_vars, region):
    """在區域 region（以 new_vars 表示的 limits 列表）上做座標變換。

    f        : 舊座標下的被積函數
    new_x    : 新座標第一個分量的表達式
    new_y    : 新座標第二個分量的表達式
    old_vars : 舊座標變數，例如 (x, y)
    new_vars : 新座標變數，例如 (r, theta)
    region   : 新座標下的積分極限，例如 [(r, 0, 1), (theta, 0, 2*pi)]
    """
    T = sp.Matrix([new_x, new_y])
    Jdet = sp.simplify(sp.det(T.jacobian(list(new_vars))))
    g = sp.simplify(f.subs(dict(zip(old_vars, [new_x, new_y]))) * sp.Abs(Jdet))
    I = sp.Integral(g, *region)          # 先建立結構，再 doit()
    return I, I.doit()

x, y = sp.symbols('x y', real=True)
r = sp.Symbol('r', positive=True, real=True)   # r > 0，讓 Abs(r) 化簡為 r
theta = sp.Symbol('theta', real=True)

# 1) 極座標變換與 Jacobian
T = sp.Matrix([r*sp.cos(theta), r*sp.sin(theta)])
J = T.jacobian([r, theta])
Jdet = sp.simplify(sp.det(J))
print('極座標 Jacobian =')
print(J)
print('det J =', Jdet)

# 2) 變數變換定理（用上面的函式實作）
# ∫∫_D f(x,y) dA = ∫∫_{T^-1(D)} f(T(r,theta)) * |det J| dr dtheta
f = x**2 + y**2
I, val = change_variables(f, r*sp.cos(theta), r*sp.sin(theta),
                          (x, y), (r, theta),
                          [(r, 0, 1), (theta, 0, 2*sp.pi)])
print('變換後的被積函數 =', I.function)
print('極座標下的積分 =', sp.simplify(val))

cart_val = sp.integrate(sp.integrate(f, (y, -sp.sqrt(1 - x**2), sp.sqrt(1 - x**2))), (x, -1, 1))
print('直角座標下的積分 =', sp.simplify(cart_val))
print('兩者相等？', sp.simplify(val - cart_val) == 0)

# 3) 圓面積：∫∫_D 1 dA（直角座標直接用對稱上下限，不需倍數因子）
area_polar = sp.integrate(sp.integrate(Jdet, (r, 0, 1)), (theta, 0, 2*sp.pi))
print('圓面積（極座標） =', area_polar)
area_cart = sp.integrate(sp.integrate(1, (y, -sp.sqrt(1 - x**2), sp.sqrt(1 - x**2))), (x, -1, 1))
print('圓面積（直角座標，對稱界限） =', area_cart)

# 對稱性的倍數要數清楚：四分之一圓盤乘 4，不是乘 2
quarter = sp.integrate(sp.integrate(1, (y, 0, sp.sqrt(1 - x**2))), (x, 0, 1))
print('四分之一圓盤 =', sp.simplify(quarter))
print('整個圓盤（乘 4）= ', sp.simplify(4*quarter))

# 4) 三重積分算球體積（極座標 det J = r^2 sin(theta)）
z = sp.Symbol('z', real=True)
phi = sp.Symbol('phi', real=True)
sph = sp.Matrix([r*sp.sin(theta)*sp.cos(phi), r*sp.sin(theta)*sp.sin(phi), r*sp.cos(theta)])
J3 = sp.simplify(sp.det(sph.jacobian([r, theta, phi])))
print('球座標 det J =', J3)
vol = sp.Integral(J3, (r, 0, 1), (theta, 0, sp.pi), (phi, 0, 2*sp.pi)).doit()
print('單位球體積 =', sp.simplify(vol))
print('對照 4/3 pi =', sp.simplify(sp.Rational(4,3)*sp.pi))

# 逐步核對 limits 順序（第一個列出者 = 內層）
s1 = sp.integrate(J3, (r, 0, 1))
s2 = sp.integrate(s1, (theta, 0, sp.pi))
s3 = sp.integrate(s2, (phi, 0, 2*sp.pi))
print('逐步 r->theta->phi =', sp.simplify(s1), sp.simplify(s2), sp.simplify(s3))

import sympy as sp

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

f = sp.exp(-(x**2 + y**2))
polar = sp.simplify(f.subs({x: r*sp.cos(theta), y: r*sp.sin(theta)}) * r)
print('極座標下的被積函數 =', polar)

val = sp.integrate(sp.integrate(polar, (r, 0, 1)), (theta, 0, 2*sp.pi))
print('∫∫_D e^{-(x^2+y^2)} dA =', sp.simplify(val))
print('= pi(1 - e^{-1}) ?', sp.simplify(val - sp.pi*(1 - sp.exp(-1))) == 0)

# 對照：直角座標下 SymPy 只能停在 erf 形式
disc = sp.integrate(sp.integrate(f, (y, -sp.sqrt(1 - x**2), sp.sqrt(1 - x**2))), (x, -1, 1))
print('直角座標結果 =', disc)
r, theta = sp.symbols('r theta', real=True)

# 整個平面上的高斯積分
I = sp.Integral(r*sp.exp(-r**2), (r, 0, sp.oo), (theta, 0, 2*sp.pi)).doit()
print('整個平面上的高斯積分 =', sp.simplify(I))

# 分開算：半徑方向與角度方向
rad = sp.integrate(r*sp.exp(-r**2), (r, 0, sp.oo))
ang = sp.integrate(1, (theta, 0, 2*sp.pi))
print('半徑方向 =', rad)
print('角度方向 =', ang)
print('乘積 =', sp.simplify(rad*ang))

# 一維高斯積分與 π 的關係（用極座標推導）
x = sp.Symbol('x', real=True)
print('一維高斯積分 =', sp.integrate(sp.exp(-x**2), (x, -sp.oo, sp.oo)))
print('二維高斯積分 / 一維的平方 =', sp.simplify(sp.pi / (sp.sqrt(sp.pi))**2))

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

# 數值逼近比較：對 ∫∫_D e^{-(x^2+y^2)} dA
def cart_riemann(n):
    xs2 = (np.arange(n) + 0.5) * (2.0/n) - 1.0
    ys2 = (np.arange(n) + 0.5) * (2.0/n) - 1.0
    X, Y = np.meshgrid(xs2, ys2)
    inside = (X**2 + Y**2) <= 1.0
    return (np.exp(-(X**2 + Y**2)) * inside).sum() * (2.0/n)**2

def polar_riemann(n):
    rs = (np.arange(n) + 0.5)/n
    ts = (np.arange(2*n) + 0.5) * (np.pi/n)
    R, T = np.meshgrid(rs, ts, indexing='ij')
    return (R*np.exp(-R**2)).sum() * (1.0/n) * (np.pi/n)

ns = [10, 20, 40, 80, 160]
for n in ns:
    print('n = %3d   cartesian = %.8f   polar = %.8f' % (n, cart_riemann(n), polar_riemann(n)))

exact = float(sp.pi*(1-sp.exp(-1)))
print('exact pi(1-e^-1) =', exact)

fig, ax = plt.subplots(figsize=(5.6, 4.2))
ax.semilogy(ns, [abs(cart_riemann(n) - exact) for n in ns], 'o-', label='Cartesian Riemann sum')
ax.semilogy(ns, [abs(polar_riemann(n) - exact) for n in ns], 's-', label='Polar Riemann sum')
ax.set_xlabel('n')
ax.set_ylabel('absolute error')
ax.set_title('Riemann sum error vs n')
ax.legend()
fig.tight_layout()
plt.show()