3.3 用 SymPy 微分：從導數定義到差分近似#

import sympy as sp

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

def deriv_by_definition(f, var, eps=0):
    """用極限定義求導數：lim_{h->0} (f(x+h) - f(x)) / h"""
    q = (f.subs(var, var + h) - f) / h
    return sp.simplify(sp.limit(q, h, eps))

# x^n
for n in [1, 2, 3, 5]:
    f = x**n
    print(f"f(x) = x^{n}:  定義法 -> {deriv_by_definition(f, x)}")

# sin
print()
print("f(x) = sin(x):  定義法 ->", deriv_by_definition(sp.sin(x), x))
print("f(x) = cos(x):  定義法 ->", deriv_by_definition(sp.cos(x), x))

# exp
print("f(x) = exp(x):  定義法 ->", deriv_by_definition(sp.exp(x), x))
print("f(x) = log(x):  定義法 ->", deriv_by_definition(sp.log(x), x))

f(x) = x^1:  定義法 -> 1
f(x) = x^2:  定義法 -> 2*x
f(x) = x^3:  定義法 -> 3*x**2
f(x) = x^5:  定義法 -> 5*x**4

f(x) = sin(x):  定義法 -> cos(x)
f(x) = cos(x):  定義法 -> -sin(x)
f(x) = exp(x):  定義法 -> exp(x)
f(x) = log(x):  定義法 -> 1/x

import sympy as sp

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

# 符號化的 x^n：n 為任意正整數
q = ((x + h)**n - x**n) / h
d = sp.limit(q, h, 0)
print("lim ((x+h)^n - x^n)/h =", d)
print("是否等於 n*x^(n-1)？", sp.simplify(d - n*x**(n - 1)) == 0)

# 對 x^n 這個函數，用 SymPy 自動展開再取極限
poly = sp.expand((x + h)**n)
print()
print("展開式（依 h 的次數）：", sp.collect(poly, h))

lim ((x+h)^n - x^n)/h = n*x**(n - 1)
是否等於 n*x^(n-1)？ True

展開式（依 h 的次數）： (h + x)**n

import sympy as sp

x = sp.Symbol('x', real=True)
Y = sp.Symbol('Y', real=True)        # 把 y 當成獨立的座標變數
u = sp.Function('u')(x)
v = sp.Function('v')(x)

# --- 乘積法則 ---
print("乘積法則： (uv)' = u'v + uv'")
lhs = sp.diff(u*v, x)
rhs = sp.diff(u, x)*v + u*sp.diff(v, x)
print("  LHS =", lhs)
print("  RHS =", rhs)
print("  相等？", sp.simplify(lhs - rhs) == 0)

# --- 鏈鎖法則 ---
print()
print("鏈鎖法則： (F(g(x)))' = F'(g(x)) g'(x)")
g = sp.Function('g')(x)
lhs2 = sp.diff(sp.sin(g), x)
rhs2 = sp.cos(g)*sp.diff(g, x)
print("  LHS =", lhs2)
print("  RHS =", rhs2)
print("  相等？", sp.simplify(lhs2 - rhs2) == 0)

# --- 隱函數微分：x^2 + y^2 = 1 ---
print()
print("隱函數微分： F(x, y) = x^2 + y^2 - 1 = 0")
F = x**2 + Y**2 - 1
Fx = sp.diff(F, x)
Fy = sp.diff(F, Y)
print("  F_x =", Fx, "   F_y =", Fy)

yp = sp.simplify(-Fx / Fy)
print("  dy/dx = -F_x/F_y =", yp)

# 與「解出 y 再微分」比對
y_explicit = sp.sqrt(1 - x**2)
print("  解出 y = sqrt(1-x^2) 再微分 =", sp.diff(y_explicit, x))
print("  代入 y 後相同？", sp.simplify(yp.subs(Y, y_explicit) - sp.diff(y_explicit, x)) == 0)

# 用 SymPy 的 Function 寫法：令 y = Function('y')(x)，但微分時把 y 當常數
print()
print("另一種寫法：用未求值的 Derivative 表達隱函數")
yy = sp.Function('yy')(x)
expr = sp.Eq(x**2 + yy**2, 1)
print("  隱函數式 :", expr)
print("  解出 y  :", sp.solve(sp.Eq(x**2 + yy**2, 1), yy))

乘積法則： (uv)' = u'v + uv'
  LHS = u(x)*Derivative(v(x), x) + v(x)*Derivative(u(x), x)
  RHS = u(x)*Derivative(v(x), x) + v(x)*Derivative(u(x), x)
  相等？ True

鏈鎖法則： (F(g(x)))' = F'(g(x)) g'(x)
  LHS = cos(g(x))*Derivative(g(x), x)
  RHS = cos(g(x))*Derivative(g(x), x)
  相等？ True

隱函數微分： F(x, y) = x^2 + y^2 - 1 = 0
  F_x = 2*x    F_y = 2*Y
  dy/dx = -F_x/F_y = -x/Y
  解出 y = sqrt(1-x^2) 再微分 = -x/sqrt(1 - x**2)
  代入 y 後相同？ True

另一種寫法：用未求值的 Derivative 表達隱函數
  隱函數式 : Eq(x**2 + yy(x)**2, 1)
  解出 y  : [-sqrt(1 - x**2), sqrt(1 - x**2)]

import numpy as np
import sympy as sp

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

# 符號真相
f = sp.sin(x)
exact = sp.diff(f, x)
x0 = float(sp.pi / 3)

exact_val = float(exact.subs(x, sp.pi / 3))
print("符號導數 f'(x) =", exact)
print("在 x0 = pi/3 ≈", x0, "處 f'(x0) =", exact_val)

def forward(f, x0, h):
    return (f(x0 + h) - f(x0)) / h

def central(f, x0, h):
    return (f(x0 + h) - f(x0 - h)) / (2 * h)

hs = [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6]
print()
print(f"{'h':>10} {'前向差商誤差':>16} {'中心差商誤差':>16}")
for h in hs:
    ef = abs(forward(np.sin, x0, h) - exact_val)
    ec = abs(central(np.sin, x0, h) - exact_val)
    print(f"{h:>10.0e} {ef:>16.3e} {ec:>16.3e}")

print()
print("誤差與 h 的比值（檢查量級）：")
for h in hs[:4]:
    ef = abs(forward(np.sin, x0, h) - exact_val)
    ec = abs(central(np.sin, x0, h) - exact_val)
    print(f"  h={h:.0e}   前向 誤差/h = {ef/h:.3e}   中心 誤差/h^2 = {ec/h**2:.3e}")

符號導數 f'(x) = cos(x)
在 x0 = pi/3 ≈ 1.0471975511965979 處 f'(x0) = 0.5

         h           前向差商誤差           中心差商誤差
     1e-01        4.410e-02        8.329e-04
     1e-02        4.338e-03        8.333e-06
     1e-03        4.331e-04        8.333e-08
     1e-04        4.330e-05        8.333e-10
     1e-05        4.330e-06        7.827e-12
     1e-06        4.330e-07        4.113e-11

誤差與 h 的比值（檢查量級）：
  h=1e-01   前向 誤差/h = 4.410e-01   中心 誤差/h^2 = 8.329e-02
  h=1e-02   前向 誤差/h = 4.338e-01   中心 誤差/h^2 = 8.333e-02
  h=1e-03   前向 誤差/h = 4.331e-01   中心 誤差/h^2 = 8.333e-02
  h=1e-04   前向 誤差/h = 4.330e-01   中心 誤差/h^2 = 8.333e-02

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

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

# 用泰勒級數從定義推出前向差商的誤差項
f = sp.sin(x)
series = sp.series(f, x, 0, 4).removeO().expand()
print("sin(x) 的三次近似 =", series)

h = sp.Symbol('h', positive=True)
q = sp.series((f.subs(x, x + h) - f) / h, h, 0, 3).removeO().expand()
print("差商 (f(x+h)-f(x))/h 的展開 =", q)
print("其中 h 的一次項 =", q.coeff(h, 1))
print("h 的二次項 =", q.coeff(h, 2))

# 數值驗證：誤差 ~ C*h
x0 = float(sp.pi / 3)
exact_val = float(sp.diff(f, x).subs(x, sp.pi / 3))

hs = np.logspace(-4, -1, 12)
err_fwd, err_ctr = [], []
for hh in hs:
    err_fwd.append(abs((np.sin(x0 + hh) - np.sin(x0)) / hh - exact_val))
    err_ctr.append(abs((np.sin(x0 + hh) - np.sin(x0 - hh)) / (2 * hh) - exact_val))

fig, ax = plt.subplots(figsize=(7, 4.5))
ax.loglog(hs, err_fwd, "o-", label="forward difference: error ~ C*h")
ax.loglog(hs, err_ctr, "s-", label="central difference: error ~ C*h^2")
ax.loglog(hs, 0.5 * hs, "k--", label="reference line C*h")
ax.loglog(hs, 0.08 * hs**2, "k:", label="reference line C*h^2")
ax.set_xlabel("h")
ax.set_ylabel("absolute error")
ax.set_title("Difference quotient error for f(x)=sin(x) at x0=pi/3")
ax.legend()
plt.tight_layout()
plt.savefig("/tmp/diff_quotient_error.png", dpi=120)
print()
print("saved /tmp/diff_quotient_error.png")
print("forward  error/h  :", [f"{e/h:.4f}" for e, h in zip(err_fwd[:4], hs[:4])])
print("central  error/h^2:", [f"{e/h**2:.4f}" for e, h in zip(err_ctr[:4], hs[:4])])

sin(x) 的三次近似 = -x**3/6 + x
差商 (f(x+h)-f(x))/h 的展開 = -h**2*cos(x)/6 - h*sin(x)/2 + cos(x)
其中 h 的一次項 = -sin(x)/2
h 的二次項 = -cos(x)/6

saved /tmp/diff_quotient_error.png
forward  error/h  : ['0.4330', '0.4330', '0.4330', '0.4331']
central  error/h^2: ['0.0833', '0.0833', '0.0833', '0.0833']

import sympy as sp

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

class Dual:
    """前向模式自動微分的一個「對數」(dual number)：值 + 一階微分"""
    def __init__(self, v, d):
        self.v, self.d = v, d
    def _pair(self, o):
        return (o.v, o.d) if isinstance(o, Dual) else (o, sp.Integer(0))
    def __add__(self, o):
        ov, od = self._pair(o)
        return Dual(self.v + ov, self.d + od)
    __radd__ = __add__
    def __sub__(self, o):
        ov, od = self._pair(o)
        return Dual(self.v - ov, self.d - od)
    def __rsub__(self, o):
        ov, od = self._pair(o)
        return Dual(ov - self.v, od - self.d)
    def __mul__(self, o):
        ov, od = self._pair(o)
        return Dual(self.v * ov, self.v * od + ov * self.d)
    __rmul__ = __mul__
    def __truediv__(self, o):
        ov, od = self._pair(o)
        return Dual(self.v / ov, (self.d * ov - self.v * od) / ov**2)
    def __rtruediv__(self, o):
        ov, od = self._pair(o)
        return Dual(ov / self.v, -ov * self.d / self.v**2)
    def __pow__(self, n):
        ov, od = self._pair(n)
        if isinstance(n, Dual):
            # 指數也是對數：用 exp(g log v) 推導
            e = n.v * sp.log(self.v)
            de = n.d * sp.log(self.v) + n.v * self.d / self.v
            return Dual(sp.exp(e), sp.exp(e) * de)
        return Dual(self.v**ov, ov * self.v**(ov - 1) * self.d)
    def __repr__(self):
        return f"Dual({self.v}, {self.d})"

def dsin(t):
    return Dual(sp.sin(t.v), sp.cos(t.v) * t.d)

def dexp(t):
    return Dual(sp.exp(t.v), sp.exp(t.v) * t.d)

def ad(f, xval):
    """回傳 (f(xval), f'(xval))"""
    t = Dual(xval, sp.Integer(1))
    r = f(t)
    return r.v, r.d

tests = [
    ("x**3 + 2*x",   lambda t: t**3 + 2*t,        x**3 + 2*x),
    ("sin(x)*exp(x)", lambda t: dsin(t) * dexp(t), sp.sin(x) * sp.exp(x)),
    ("x**x",          lambda t: t**t,             x**x),
    ("1/x",           lambda t: 1 / t,            1 / x),
]

for name, adf, expr in tests:
    val, der = ad(adf, sp.Integer(2))
    sym = sp.diff(expr, x).subs(x, 2)
    print(f"{name:>16}:  AD 值 = {val},  AD 導數 = {der},  SymPy 導數 = {sym},  相等？ {sp.simplify(der - sym) == 0}")

      x**3 + 2*x:  AD 值 = 12,  AD 導數 = 14,  SymPy 導數 = 14,  相等？ True
   sin(x)*exp(x):  AD 值 = exp(2)*sin(2),  AD 導數 = exp(2)*cos(2) + exp(2)*sin(2),  SymPy 導數 = exp(2)*cos(2) + exp(2)*sin(2),  相等？ True
            x**x:  AD 值 = 4,  AD 導數 = 4*log(2) + 4,  SymPy 導數 = 4*log(2) + 4,  相等？ True
             1/x:  AD 值 = 1/2,  AD 導數 = -1/4,  SymPy 導數 = -1/4,  相等？ True

import numpy as np
import sympy as sp

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

def deriv_by_definition(f, var):
    return sp.simplify(sp.limit((f.subs(var, var + h) - f) / h, h, 0))

# 一個略帶挑戰的函數：x*exp(-x^2) 的導數
f = x * sp.exp(-x**2)
print("f(x) = x*exp(-x^2)")
print("  定義法  ->", deriv_by_definition(f, x))
print("  sp.diff ->", sp.simplify(sp.diff(f, x)))
print("  一致？", sp.simplify(deriv_by_definition(f, x) - sp.diff(f, x)) == 0)

# 數值檢查
x0 = 1.3
fd = sp.lambdify(x, sp.diff(f, x), "numpy")
exact = float(fd(x0))
print()
print("x0 =", x0, " 精確導數值 =", exact)

for hh in [1e-2, 1e-4, 1e-6]:
    num = (float(f.subs(x, x0 + hh)) - float(f.subs(x, x0))) / hh
    print(f"  h = {hh:.0e}  差商 = {num:.12f}  誤差 = {abs(num - exact):.3e}")

# 極值點：費馬原理的實用版
print()
print("f'(x) = 0 的解（極值點）= ", sp.solve(sp.Eq(sp.diff(f, x), 0), x))
for c in sp.solve(sp.Eq(sp.diff(f, x), 0), x):
    print("  x =", c, " f(x) =", sp.simplify(f.subs(x, c)))

f(x) = x*exp(-x^2)
  定義法  -> (1 - 2*x**2)*exp(-x**2)
  sp.diff -> (1 - 2*x**2)*exp(-x**2)
  一致？ True

x0 = 1.3  精確導數值 = -0.43915646710331446
  h = 1e-02  差商 = -0.438209220029  誤差 = 9.472e-04
  h = 1e-04  差商 = -0.439147348237  誤差 = 9.119e-06
  h = 1e-06  差商 = -0.439156375931  誤差 = 9.117e-08

f'(x) = 0 的解（極值點）=  [-sqrt(2)/2, sqrt(2)/2]
  x = -sqrt(2)/2  f(x) = -sqrt(2)*exp(-1/2)/2
  x = sqrt(2)/2  f(x) = sqrt(2)*exp(-1/2)/2