10.3 用 SymPy 做複變數：留數、展開與圍線積分#

import sympy as sp

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

print("=== 1. SymPy 的複數常數 ===")
print("  sp.I =", sp.I, " ; sp.I**2 =", sp.I**2)
print("  sp.re(3+4*I) =", sp.re(3 + 4 * sp.I))
print("  sp.im(3+4*I) =", sp.im(3 + 4 * sp.I))
print("  sp.Abs(3+4*I) =", sp.Abs(3 + 4 * sp.I))
print("  sp.arg(3+4*I) =", sp.arg(3 + 4 * sp.I))
print("  sp.conjugate(3+4*I) =", sp.conjugate(3 + 4 * sp.I))
print("  (3+4*I)/(3-4*I) =", sp.simplify((3 + 4 * sp.I) / (3 - 4 * sp.I)))
print()

print("=== 2. exp / sin / cos 在複數上的行為 ===")
w = sp.exp(sp.I * sp.pi / 3)
print("  exp(i*pi/3) =", w)
print("  expand_complex =", sp.expand_complex(w))
print("  |exp(i*pi/3)| =", sp.Abs(w), " ; arg =", sp.arg(w))
print("  sin(pi/6 + i*pi/4) =", sp.expand_complex(sp.sin(sp.pi / 6 + sp.I * sp.pi / 4)))
print("  cos(pi/6 + i*pi/4) =", sp.expand_complex(sp.cos(sp.pi / 6 + sp.I * sp.pi / 4)))
print("  exp(1+2i) =", sp.expand_complex(sp.exp(1 + 2 * sp.I)))
print("  exp(1+2i) 的絕對值 =", sp.Abs(sp.exp(1 + 2 * sp.I)), " (= e^1 =", sp.E, ")")
print()

print("=== 3. expand_complex 與 re/im ===")
expr = (1 + sp.I) ** 5
print("  (1+i)^5 =", expr)
print("  expand_complex((1+i)^5) =", sp.expand_complex(expr))
print("  re =", sp.re(expr), " ; im =", sp.im(expr))
print()
f = sp.Function('f')
expr2 = sp.exp(sp.I * x) * sp.cos(y)
print("  f = exp(i*x)*cos(y)")
print("  re(f) =", sp.simplify(sp.re(expr2)))
print("  im(f) =", sp.simplify(sp.im(expr2)))
print("  expand_complex(f) =", sp.expand_complex(expr2))
print()

print("=== 4. sp.re / sp.im 對實符號與假設的處理 ===")
print("  re(sqrt(x^2+y^2)) =", sp.re(sp.sqrt(x**2 + y**2)))
print("  im(log(z)) (z 為複數符號) =", sp.im(sp.log(z)))
print("  simplify(expand_complex(log(z))) =", sp.simplify(sp.expand_complex(sp.log(z))))

=== 1. SymPy 的複數常數 ===
  sp.I = I  ; sp.I**2 = -1
  sp.re(3+4*I) = 3
  sp.im(3+4*I) = 4
  sp.Abs(3+4*I) = 5
  sp.arg(3+4*I) = atan(4/3)
  sp.conjugate(3+4*I) = 3 - 4*I
  (3+4*I)/(3-4*I) = -7/25 + 24*I/25

=== 2. exp / sin / cos 在複數上的行為 ===
  exp(i*pi/3) = exp(I*pi/3)
  expand_complex = 1/2 + sqrt(3)*I/2
  |exp(i*pi/3)| = 1  ; arg = pi/3
  sin(pi/6 + i*pi/4) = cosh(pi/4)/2 + sqrt(3)*I*sinh(pi/4)/2
  cos(pi/6 + i*pi/4) = sqrt(3)*cosh(pi/4)/2 - I*sinh(pi/4)/2
  exp(1+2i) = E*cos(2) + E*I*sin(2)
  exp(1+2i) 的絕對值 = E  (= e^1 = E )

=== 3. expand_complex 與 re/im ===
  (1+i)^5 = (1 + I)**5
  expand_complex((1+i)^5) = -4 - 4*I
  re = -4  ; im = -4

  f = exp(i*x)*cos(y)
  re(f) = cos(x)*cos(y)
  im(f) = sin(x)*cos(y)
  expand_complex(f) = I*sin(x)*cos(y) + cos(x)*cos(y)

=== 4. sp.re / sp.im 對實符號與假設的處理 ===
  re(sqrt(x^2+y^2)) = sqrt(x**2 + y**2)
  im(log(z)) (z 為複數符號) = arg(z)
  simplify(expand_complex(log(z))) = log(re(z)**2 + im(z)**2)/2 + I*arg(re(z) + I*im(z))

import sympy as sp

z = sp.Symbol('z')
w = sp.Symbol('w')

print("=== 1. 在複數點做 Taylor 展開 ===")
a = 1 + sp.I
f = sp.exp(z) / (z**2 + 1)
s = sp.series(f, z, a, 3)
print("  exp(z)/(z^2+1) 在 a = 1+i 的三次展開：")
print("  ", sp.simplify(s))
print("  截斷誤差 =", sp.simplify(s.removeO().subs(z, a) - f.subs(z, a)))
print()

print("=== 2. Laurent 展開（含負冪次）===")
print("  1/(z^2-1) 在 z=1 的 Laurent 展開 =", sp.series(1 / (z**2 - 1), z, 1, 2))
print("  1/(z-1) 在 z=1 的 Laurent 展開   =", sp.series(1 / (z - 1), z, 1, 1))
print("  z/(z-1)^2 在 z=1 的 Laurent 展開 =", sp.series(z / (z - 1)**2, z, 1, 2))
print()

print("=== 3. 留數 = Laurent 展開中 (z-a)^{-1} 的係數 ===")
cases = [(1 / (z**2 + 1), sp.I), (sp.exp(z) / z**2, 0), (sp.sin(z) / z**3, 0),
         (z / (z - 1)**2, 1)]
for expr, a0 in cases:
    s2 = sp.series(expr, z, a0, 1)
    # 令 z = a0 + w，把展開改寫成 w 的冪級數，再取出 w^{-1} 的係數
    shifted = sp.expand(s2.removeO().subs(z, a0 + w))
    coef = sp.simplify(sp.expand(shifted).coeff(w, -1))
    res = sp.residue(expr, z, a0)
    print("  f(z) = %-22s 在 z=%-4s 的展開 = %s" % (expr, a0, sp.series(expr, z, a0, 1)))
    print("      展開中的 (z-a)^(-1) 係數 = %s ； sp.residue = %s ； 相同？ %s"
          % (coef, res, sp.simplify(coef - res) == 0))
print()

print("=== 4. 收斂半徑由最近的奇點決定 ===")
h = 1 / ((z - 1) * (z + 2))
print("  f(z) = 1/((z-1)(z+2))，部分分式 =", sp.apart(h, z, extension=sp.I))
print("  在 0 的 Taylor 展開 =", sp.series(h, z, 0, 3))
n = sp.symbols('n', integer=True, nonnegative=True)
# 展開式：-(1/3) * sum z^n  -  (1/6) * sum (-z/2)^n
def partial_sum(zval, N):
    s1 = sp.summation(zval**n, (n, 0, N))
    s2 = sp.summation((-zval / 2)**n, (n, 0, N))
    return sp.simplify(-sp.Rational(1, 3) * s1 - sp.Rational(1, 6) * s2)

print()
for zval, tag in [(sp.Rational(1, 2), "收斂"), (sp.Rational(3, 2), "發散")]:
    print("  z = %s（%s）：" % (zval, tag))
    for N in [5, 100]:
        print("      N=%-4d 部分和 = %s" % (N, sp.N(partial_sum(zval, N), 12)))
print("  真值 f(0.5) =", sp.N(h.subs(z, sp.Rational(1, 2)), 12))

=== 1. 在複數點做 Taylor 展開 ===
  exp(z)/(z^2+1) 在 a = 1+i 的三次展開：
   (152 - 164*I)*exp(1 + I)/250 + z*(-34 + 63*I)*exp(1 + I)/125 + z**2*(3 - 46*I)*exp(1 + I)/250 + O((z - 1 - I)**3, (z, 1 + I))
  截斷誤差 = 0

=== 2. Laurent 展開（含負冪次）===
  1/(z^2-1) 在 z=1 的 Laurent 展開 = 1/(2*(z - 1)) - 3/8 + z/8 + O((z - 1)**2, (z, 1))
  1/(z-1) 在 z=1 的 Laurent 展開   = 1/(z - 1)
  z/(z-1)^2 在 z=1 的 Laurent 展開 = (z - 1)**(-2) + 1/(z - 1) + O((z - 1)**2, (z, 1))

=== 3. 留數 = Laurent 展開中 (z-a)^{-1} 的係數 ===
  f(z) = 1/(z**2 + 1)           在 z=I    的展開 = -I/(2*(z - I)) + 1/4 + O(z - I, (z, I))
      展開中的 (z-a)^(-1) 係數 = -I/2 ； sp.residue = -I/2 ； 相同？ True
  f(z) = exp(z)/z**2            在 z=0    的展開 = z**(-2) + 1/z + 1/2 + O(z)
      展開中的 (z-a)^(-1) 係數 = 1 ； sp.residue = 1 ； 相同？ True
  f(z) = sin(z)/z**3            在 z=0    的展開 = z**(-2) - 1/6 + O(z)
      展開中的 (z-a)^(-1) 係數 = 0 ； sp.residue = 0 ； 相同？ True
  f(z) = z/(z - 1)**2           在 z=1    的展開 = (z - 1)**(-2) + 1/(z - 1) + O(z - 1, (z, 1))
      展開中的 (z-a)^(-1) 係數 = 1 ； sp.residue = 1 ； 相同？ True

=== 4. 收斂半徑由最近的奇點決定 ===
  f(z) = 1/((z-1)(z+2))，部分分式 = -1/(3*(z + 2)) + 1/(3*(z - 1))
  在 0 的 Taylor 展開 = -1/2 - z/4 - 3*z**2/8 + O(z**3)

  z = 1/2（收斂）：
      N=5    部分和 = -0.789550781250
      N=100  部分和 = -0.800000000000
  z = 3/2（發散）：
      N=5    部分和 = -7.00537109375
      N=100  部分和 = -4.06561177535e+17
  真值 f(0.5) = -0.800000000000

import sympy as sp

z = sp.Symbol('z')
x = sp.Symbol('x', real=True)
w = sp.Symbol('w')

def upper_poles(g):
    """找出 g(x) 作為複變數函數時，上半平面的極點。"""
    gs = g.subs(x, z)
    num, den = sp.fraction(sp.together(gs))
    roots = sp.solve(sp.Eq(den, 0), z)
    out = []
    for r in roots:
        im = float(sp.im(sp.N(r)))
        if im > 1e-12 and sp.simplify(num.subs(z, r)) != 0:
            out.append(r)
    return out

def residue_integral(g):
    """用留數法計算 ∫_{-∞}^{∞} g(x) dx = 2πi Σ_{上半平面} Res。"""
    poles = upper_poles(g)
    total = sum(sp.residue(g.subs(x, z), z, p) for p in poles)
    return sp.simplify(sp.expand_complex(2 * sp.pi * sp.I * total)), poles

def residue_at_infinity(f, zsym=z):
    """Res(∞, f) = -Res_{w=0} f(1/w)/w^2（SymPy 尚未實作 sp.residue(f, z, sp.oo)）。"""
    return -sp.residue(sp.simplify(f.subs(zsym, 1 / w) / w**2), w, 0)

print("=== 1. sp.residue：各種極點的留數 ===")
tests = [(1 / z, 0), (1 / z**2, 0), (sp.sin(z) / z**3, 0),
         (sp.exp(z) / z**2, 0), (z / (z - 1)**3, 1), (1 / (z**2 + 1), sp.I)]
for f, a0 in tests:
    print("  f(z) = %-16s 在 z=%-4s 的留數 = %s" % (f, a0, sp.residue(f, z, a0)))
print()

print("=== 2. 留數法算 ∫_{-∞}^{∞} dx（g(z) 在上半平面有界時可用）===")
integrands = [
    ("1/(x**2+1)",          1 / (x**2 + 1)),
    ("1/(x**2+2)",          1 / (x**2 + 2)),
    ("1/(x**2-2*x+5)",      1 / (x**2 - 2 * x + 5)),
    ("1/(x**4+1)",          1 / (x**4 + 1)),
    ("1/(x**2+1)**2",       1 / (x**2 + 1)**2),
    ("x/(x**2+1)**2",       x / (x**2 + 1)**2),
]
for name, g in integrands:
    val, poles = residue_integral(g)
    direct = sp.simplify(sp.integrate(g, (x, -sp.oo, sp.oo)))
    print("  ∫ %-20s 留數法 = %-20s 直接積分 = %-20s 相等？ %-5s 極點 = %s"
          % (name + " dx", val, direct, sp.simplify(val - direct) == 0, poles))
print()

print("=== 3. 陷阱：把 cos 直接送進留數法 ===")
g = sp.cos(x) / (x**2 + 1)
naive, poles = residue_integral(g)
direct = sp.integrate(g, (x, -sp.oo, sp.oo))
print("  留數法（直接用 cos）= %s ； 正確答案 = %s ； 相等？ %s"
      % (naive, sp.simplify(direct), sp.simplify(naive - direct) == 0))
print("  原因：cos(z) = (e^{iz}+e^{-iz})/2，e^{-iz} 在上半平面是爆炸的，")
print("        大圓弧上的積分不會趨近 0，所以不能直接套上半平面閉合。")
print()
print("  正確做法：用 e^{iz}（上半平面衰減），再取實部")
F = sp.exp(sp.I * z) / (z**2 + 1)
res = sp.residue(F, z, sp.I)
print("    Res_{z=i} e^{iz}/(z^2+1) =", sp.simplify(res))
print("    2πi * Res =", sp.simplify(sp.expand_complex(2 * sp.pi * sp.I * res)))
print("    其實部即 ∫ cos(x)/(x^2+1) dx =",
      sp.simplify(sp.re(2 * sp.pi * sp.I * res)), " ； 直接積分 =", sp.simplify(direct))
print()

print("=== 4. 留數定理：∮ f dz = 2πi Σ Res ===")
f = 1 / (z**2 + 1)
print("  ∮_{|z|=2} dz/(z^2+1) = 2πi (Res(i) + Res(-i)) =",
      sp.simplify(2 * sp.pi * sp.I * (sp.residue(f, z, sp.I) + sp.residue(f, z, -sp.I))))
print("  ∮_{|z-i|=1/2} dz/(z^2+1) = 2πi Res(i) =",
      sp.simplify(2 * sp.pi * sp.I * sp.residue(f, z, sp.I)))
print()

print("=== 5. 無窮遠點的留數 ===")
for name, g in [("1/(z**2+1)", 1 / (z**2 + 1)), ("1/z", 1 / z), ("1/(z-1)", 1 / (z - 1))]:
    finite = sp.solve(sp.denom(sp.together(g)), z)
    s = sum(sp.residue(g, z, p) for p in finite)
    print("  %-14s 有限留數之和 = %-8s Res(∞) = %-8s 相加為 0？ %s"
          % (name, sp.simplify(s), sp.simplify(residue_at_infinity(g)),
             sp.simplify(s + residue_at_infinity(g)) == 0))

=== 1. sp.residue：各種極點的留數 ===
  f(z) = 1/z              在 z=0    的留數 = 1
  f(z) = z**(-2)          在 z=0    的留數 = 0
  f(z) = sin(z)/z**3      在 z=0    的留數 = 0
  f(z) = exp(z)/z**2      在 z=0    的留數 = 1
  f(z) = z/(z - 1)**3     在 z=1    的留數 = 0
  f(z) = 1/(z**2 + 1)     在 z=I    的留數 = -I/2

=== 2. 留數法算 ∫_{-∞}^{∞} dx（g(z) 在上半平面有界時可用）===
  ∫ 1/(x**2+1) dx        留數法 = pi                   直接積分 = pi                   相等？ True  極點 = [I]
  ∫ 1/(x**2+2) dx        留數法 = sqrt(2)*pi/2         直接積分 = sqrt(2)*pi/2         相等？ True  極點 = [sqrt(2)*I]
  ∫ 1/(x**2-2*x+5) dx    留數法 = pi/2                 直接積分 = pi/2                 相等？ True  極點 = [1 + 2*I]
  ∫ 1/(x**4+1) dx        留數法 = sqrt(2)*pi/2         直接積分 = sqrt(2)*pi/2         相等？ True  極點 = [-sqrt(2)/2 + sqrt(2)*I/2, sqrt(2)/2 + sqrt(2)*I/2]
  ∫ 1/(x**2+1)**2 dx     留數法 = pi/2                 直接積分 = pi/2                 相等？ True  極點 = [I]
  ∫ x/(x**2+1)**2 dx     留數法 = 0                    直接積分 = 0                    相等？ True  極點 = [I]

=== 3. 陷阱：把 cos 直接送進留數法 ===
  留數法（直接用 cos）= pi*cosh(1) ； 正確答案 = pi*exp(-1) ； 相等？ False
  原因：cos(z) = (e^{iz}+e^{-iz})/2，e^{-iz} 在上半平面是爆炸的，
        大圓弧上的積分不會趨近 0，所以不能直接套上半平面閉合。

  正確做法：用 e^{iz}（上半平面衰減），再取實部
    Res_{z=i} e^{iz}/(z^2+1) = -I*exp(-1)/2
    2πi * Res = pi*exp(-1)
    其實部即 ∫ cos(x)/(x^2+1) dx = pi*exp(-1)  ； 直接積分 = pi*exp(-1)

=== 4. 留數定理：∮ f dz = 2πi Σ Res ===
  ∮_{|z|=2} dz/(z^2+1) = 2πi (Res(i) + Res(-i)) = 0
  ∮_{|z-i|=1/2} dz/(z^2+1) = 2πi Res(i) = pi

=== 5. 無窮遠點的留數 ===
  1/(z**2+1)     有限留數之和 = 0        Res(∞) = 0        相加為 0？ True
  1/z            有限留數之和 = 1        Res(∞) = -1       相加為 0？ True
  1/(z-1)        有限留數之和 = 1        Res(∞) = -1       相加為 0？ True

import sympy as sp

z = sp.Symbol('z')
t = sp.Symbol('t', real=True)

print("=== 1. 用路徑參數化計算複積分 ===")
# (a) 線段 z(t) = t, 0 <= t <= 1：從 0 到 1
I1 = sp.integrate(t**2, (t, 0, 1))
print("  (a) ∫ z^2 dz 沿直線段 z = t （0 -> 1）     =", I1, " ； 理論 1/3 ； 相等？", I1 == sp.Rational(1, 3))

# (b) 直線段 z(t) = t(1+i), 0 <= t <= 1：從 0 到 1+i
zt = t * (1 + sp.I)
dz = sp.diff(zt, t)
I2 = sp.simplify(sp.integrate(zt**2 * dz, (t, 0, 1)))
print("  (b) ∫ z^2 dz 沿直線段 z = t(1+i)（0 -> 1+i）=", I2, " ； 理論 (1+i)^3/3 ； 相等？",
      sp.simplify(I2 - (1 + sp.I)**3 / 3) == 0)

# (c) 彎曲路徑：起訖點同為 0 -> 1+i
zt2 = t * (1 + sp.I) + sp.I * t * (1 - t)
dz2 = sp.diff(zt2, t)
I3 = sp.simplify(sp.integrate(zt2**2 * dz2, (t, 0, 1)))
print("  (c) ∫ z^2 dz 沿彎曲路徑 z = t(1+i)+i t(1-t)  =", I3, " ； 與 (b) 相等？",
      sp.simplify(I3 - I2) == 0)
print("      （全純函數的路徑積分只取決於起訖點，與路徑無關 —— 科西積分定理的直接後果）")
print()

print("=== 2. 圓周上的參數化積分：∮ dz/z = 2πi ===")
R = sp.Integer(1)
a = sp.Integer(0)
zt = a + R * sp.exp(sp.I * t)
integrand = sp.simplify((1 / zt) * sp.diff(zt, t))
print("  z = e^{it}，dz/dt = i e^{it}，被積函數 =", integrand)
I4 = sp.simplify(sp.integrate(integrand, (t, 0, 2 * sp.pi)))
print("  ∮_{|z|=1} dz/z =", I4, " ； 等於 2πi？ ", I4 == 2 * sp.pi * sp.I)
print()

print("=== 3. 用 sp.Integral 保留未求值形式，再手動 .doit() ===")
expr = sp.Integral(sp.exp(z), (z, 0, 1))
print("  保留形式 =", expr)
print("  .doit() =", expr.doit())
print()

print("=== 4. 路徑無關性的另一個證明：對數 ===")
print("  ∫_{1}^{-1} dz/z 沿上半平面圓弧（z=e^{it}, t: 0->π）:")
zt3 = sp.exp(sp.I * t)
print("     ", sp.simplify(sp.integrate((1 / zt3) * sp.diff(zt3, t), (t, 0, sp.pi))))
print("  ∫_{1}^{-1} dz/z 沿下半平面圓弧（z=e^{-it}, t: 0->π）:")
zt4 = sp.exp(-sp.I * t)
print("     ", sp.simplify(sp.integrate((1 / zt4) * sp.diff(zt4, t), (t, 0, sp.pi))))
print("  => 相差 2πi，正是 log 的多值性（分支切割）。")

=== 1. 用路徑參數化計算複積分 ===
  (a) ∫ z^2 dz 沿直線段 z = t （0 -> 1）     = 1/3  ； 理論 1/3 ； 相等？ True
  (b) ∫ z^2 dz 沿直線段 z = t(1+i)（0 -> 1+i）= -2/3 + 2*I/3  ； 理論 (1+i)^3/3 ； 相等？ True
  (c) ∫ z^2 dz 沿彎曲路徑 z = t(1+i)+i t(1-t)  = -2/3 + 2*I/3  ； 與 (b) 相等？ True
      （全純函數的路徑積分只取決於起訖點，與路徑無關 —— 科西積分定理的直接後果）

=== 2. 圓周上的參數化積分：∮ dz/z = 2πi ===
  z = e^{it}，dz/dt = i e^{it}，被積函數 = I
  ∮_{|z|=1} dz/z = 2*I*pi  ； 等於 2πi？  True

=== 3. 用 sp.Integral 保留未求值形式，再手動 .doit() ===
  保留形式 = Integral(exp(z), (z, 0, 1))
  .doit() = -1 + E

=== 4. 路徑無關性的另一個證明：對數 ===
  ∫_{1}^{-1} dz/z 沿上半平面圓弧（z=e^{it}, t: 0->π）:
      I*pi
  ∫_{1}^{-1} dz/z 沿下半平面圓弧（z=e^{-it}, t: 0->π）:
      -I*pi
  => 相差 2πi，正是 log 的多值性（分支切割）。

import sympy as sp

z = sp.Symbol('z')
w = sp.Symbol('w')

print("=== 1. Cauchy 公式：f^(n)(a) = n!/(2πi) ∮ f(z)/(z-a)^{n+1} dz ===")
a = sp.Rational(3, 10) + sp.I * sp.Rational(1, 5)
t = sp.Symbol('t', real=True)
R = sp.Rational(1, 2)
zt = a + R * sp.exp(sp.I * t)
dz = sp.diff(zt, t)

fpoly = z**3 + 2 * z
for n in range(0, 4):
    integ = sp.expand(fpoly.subs(z, zt) / (zt - a) ** (n + 1) * dz)
    val = sp.simplify(sp.expand_complex(sp.integrate(integ, (t, 0, 2 * sp.pi))))
    target = 2 * sp.pi * sp.I * sp.diff(fpoly, z, n).subs(z, a) / sp.factorial(n)
    print("  n = %d : ∮ f(z)/(z-a)^{n+1} dz = %-28s 與 2πi f^(n)(a)/n! 之差 = %s"
          % (n, val, sp.simplify(val - sp.expand_complex(target))))
print()

print("=== 2. 同一件事的另一種看法：把 f 展開成冪級數再逐項積分 ===")
for k in range(0, 4):
    ck = sp.simplify(sp.expand_complex(sp.diff(fpoly, z, k).subs(z, a) / sp.factorial(k)))
    contrib = "2πi" if k == 0 else "0"
    print("      k = %d : 係數 c_%d = %-24s ； ∮ (z-a)^{k-1} dz = %s" % (k, k, ck, contrib))
print("  => 只有 k = 0 那一項貢獻 2πi，其餘貢獻 0，這正是 f(a)。")
print("  對 n > 0 同理：只有 k = n 貢獻，得到 2πi f^(n)(a)/n!。")
print()

print("=== 3. 用留數定理重寫 Cauchy 公式 ===")
g = fpoly / (z - a) ** 4
print("  g(z) = f(z)/(z-a)^4 =", g)
print("  Res(g, a) =", sp.simplify(sp.residue(g, z, a)))
print("  f^(3)(a)/3! =", sp.simplify(sp.diff(fpoly, z, 3).subs(z, a) / sp.factorial(3)))
print("  兩者相等？",
      sp.simplify(sp.residue(g, z, a) - sp.diff(fpoly, z, 3).subs(z, a) / sp.factorial(3)) == 0)
print()

print("=== 4. 一個實用的副產品：求 Taylor 係數 ===")
c2 = sp.simplify(sp.expand_complex(sp.diff(fpoly, z, 2).subs(z, a) / sp.factorial(2)))
print("  f^(2)(a)/2! =", c2)
series_part = sp.expand(sp.series(fpoly, z, a, 3).removeO().subs(z, a + w))
print("  用 series 展開後取 w^2 係數 =", sp.simplify(series_part.coeff(w, 2)))
print("  兩者相等？", sp.simplify(c2 - series_part.coeff(w, 2)) == 0)

=== 1. Cauchy 公式：f^(n)(a) = n!/(2πi) ∮ f(z)/(z-a)^{n+1} dz ===
  n = 0 : ∮ f(z)/(z-a)^{n+1} dz = pi*(-446 + 591*I)/500        與 2πi f^(n)(a)/n! 之差 = 0
  n = 1 : ∮ f(z)/(z-a)^{n+1} dz = pi*(-36 + 215*I)/50          與 2πi f^(n)(a)/n! 之差 = 0
  n = 2 : ∮ f(z)/(z-a)^{n+1} dz = pi*(-6 + 9*I)/5              與 2πi f^(n)(a)/n! 之差 = 0
  n = 3 : ∮ f(z)/(z-a)^{n+1} dz = 2*I*pi                       與 2πi f^(n)(a)/n! 之差 = 0

=== 2. 同一件事的另一種看法：把 f 展開成冪級數再逐項積分 ===
      k = 0 : 係數 c_0 = 591/1000 + 223*I/500     ； ∮ (z-a)^{k-1} dz = 2πi
      k = 1 : 係數 c_1 = 43/20 + 9*I/25           ； ∮ (z-a)^{k-1} dz = 0
      k = 2 : 係數 c_2 = 9/10 + 3*I/5             ； ∮ (z-a)^{k-1} dz = 0
      k = 3 : 係數 c_3 = 1                        ； ∮ (z-a)^{k-1} dz = 0
  => 只有 k = 0 那一項貢獻 2πi，其餘貢獻 0，這正是 f(a)。
  對 n > 0 同理：只有 k = n 貢獻，得到 2πi f^(n)(a)/n!。

=== 3. 用留數定理重寫 Cauchy 公式 ===
  g(z) = f(z)/(z-a)^4 = (z**3 + 2*z)/(z - 3/10 - I/5)**4
  Res(g, a) = 1
  f^(3)(a)/3! = 1
  兩者相等？ True

=== 4. 一個實用的副產品：求 Taylor 係數 ===
  f^(2)(a)/2! = 9/10 + 3*I/5
  用 series 展開後取 w^2 係數 = 9/10 + 3*I/5
  兩者相等？ True