13.3 用 SymPy 做符號最佳化與自動微分#

import sympy as sp

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

# ============================================================
# 1. 用 sp.diff + sp.solve 求臨界點
# ============================================================
print("=== 1. 用 sp.diff + sp.solve 求臨界點 ===")
f = x**4 - 2*x**2 + 1
df = sp.diff(f, x)
print("f  =", f)
print("f' =", df)
crit = sp.solve(sp.Eq(df, 0), x)
print("f' = 0 的解 =", crit)

# ============================================================
# 2. 二階判別
# ============================================================
print()
print("=== 2. 二階判別 ===")
d2f = sp.diff(f, x, 2)
print("f'' =", d2f)
for c in crit:
    v2 = sp.simplify(d2f.subs(x, c))
    fv = sp.simplify(f.subs(x, c))
    kind = "極小" if v2 > 0 else ("極大" if v2 < 0 else "不確定")
    print(f"   x = {c}:  f'' = {v2}   f = {fv}   -> {kind}")

# ============================================================
# 3. sp.minimum / sp.maximum
# ============================================================
print()
print("=== 3. sp.minimum / sp.maximum ===")
print("minimum(x**2+1, x) =", sp.minimum(x**2 + 1, x))
print("maximum(-x**2, x)  =", sp.maximum(-x**2, x))
print("minimum(sin(x), x)  =", sp.minimum(sp.sin(x), x))
print("maximum(x**3 - 3*x, x) =", sp.maximum(x**3 - 3*x, x))

# ============================================================
# 4. Hessian 判別正定性（Sylvester 準則）
# ============================================================
print()
print("=== 4. Hessian 判別正定性 ===")
y = sp.symbols('y', real=True)
g = 4*x**2 + 2*y**2 + 4*x*y - 4*x - 2*y
H = sp.hessian(g, (x, y))
print("g  =", g)
print("grad =", (sp.diff(g, x), sp.diff(g, y)))
crit2 = sp.solve([sp.diff(g, x), sp.diff(g, y)], (x, y), dict=True)
print("臨界點 =", crit2)
Hc = H.subs(crit2[0])
print("Hessian at crit =", Hc.tolist())
print("順序主元（Sylvester）：", [sp.det(Hc[:k, :k]) for k in range(1, 3)])
print("g(臨界點) =", sp.simplify(g.subs(crit2[0])))

=== 1. 用 sp.diff + sp.solve 求臨界點 ===
f  = x**4 - 2*x**2 + 1
f' = 4*x**3 - 4*x
f' = 0 的解 = [-1, 0, 1]

=== 2. 二階判別 ===
f'' = 4*(3*x**2 - 1)
   x = -1:  f'' = 8   f = 0   -> 極小
   x = 0:  f'' = -4   f = 1   -> 極大
   x = 1:  f'' = 8   f = 0   -> 極小

=== 3. sp.minimum / sp.maximum ===
minimum(x**2+1, x) = 1
maximum(-x**2, x)  = 0
minimum(sin(x), x)  = -1
maximum(x**3 - 3*x, x) = oo

=== 4. Hessian 判別正定性 ===
g  = 4*x**2 + 4*x*y - 4*x + 2*y**2 - 2*y
grad = (8*x + 4*y - 4, 4*x + 4*y - 2)
臨界點 = [{x: 1/2, y: 0}]
Hessian at crit = [[8, 4], [4, 4]]
順序主元（Sylvester）： [8, 16]
g(臨界點) = -1

import sympy as sp

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

# 簡單的多項式：SymPy 能給出封閉答案
print("minimum(x**2+1, x)      =", sp.minimum(x**2 + 1, x))
print("minimum((x-2)**2, x)    =", sp.minimum((x - 2)**2, x))
print("minimum(sin(x)*cos(x),x)=", sp.minimum(sp.sin(x) * sp.cos(x), x))

# 超越方程：SymPy 直接承認失敗，而不是編一個數字
try:
    print("minimum(sin(x)+x**3, x)=", sp.minimum(sp.sin(x) + x**3, x))
except NotImplementedError as e:
    print("minimum(sin(x)+x**3, x) 失敗：", e)

minimum(x**2+1, x)      = 1
minimum((x-2)**2, x)    = 0
minimum(sin(x)*cos(x),x)= -1/2
minimum(sin(x)+x**3, x) 失敗： Unable to find critical points for x**3 + sin(x)

import sympy as sp

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

# ============================================================
# 案例 1：最大化 x*y，約束 x^2 + y^2 = 1
# ============================================================
print("=== 案例 1：max x*y  s.t.  x^2 + y^2 = 1 ===")
f = x*y
g = x**2 + y**2 - 1
L = f - lam*g
print("目標 f =", f)
print("約束 g =", g)
print("拉格朗日函數 L =", L)

eqs = [sp.diff(L, v) for v in (x, y)] + [g]
print("特點方程組：")
for e in eqs:
    print("   ", e)

sols = sp.solve(eqs, (x, y, lam), dict=True)
print()
print("解：")
for s in sols:
    print("   ", {k: sp.simplify(v) for k, v in s.items()}, "   f =", sp.simplify(f.subs(s)))

# ============================================================
# 案例 2：最大化 x + y，約束 x^2 + y^2 = 1
# ============================================================
print()
print("=== 案例 2：max x + y  s.t.  x^2 + y^2 = 1 ===")
f2 = x + y
L2 = f2 - lam*g
eqs2 = [sp.diff(L2, v) for v in (x, y)] + [g]
for s in sp.solve(eqs2, (x, y, lam), dict=True):
    print("   ", {k: sp.simplify(v) for k, v in s.items()}, "  f =", sp.simplify(f2.subs(s)))

# ============================================================
# 案例 3：帶不等式約束（用 KKT 的形式）
# ============================================================
print()
print("=== 案例 3：min x^2 + y^2  s.t.  x + y >= 2 ===")
f3 = x**2 + y**2
g3 = x + y - 2              # 約束 g3 >= 0
L3 = f3 - lam*g3            # 對 >= 約束取 - 號，並要求 lam >= 0（KKT 形式）
eqs3 = [sp.diff(L3, v) for v in (x, y)] + [g3]
for s in sp.solve(eqs3, (x, y, lam), dict=True):
    print("   ", {k: sp.simplify(v) for k, v in s.items()},
          " f =", sp.simplify(f3.subs(s)))

=== 案例 1：max x*y  s.t.  x^2 + y^2 = 1 ===
目標 f = x*y
約束 g = x**2 + y**2 - 1
拉格朗日函數 L = -lam*(x**2 + y**2 - 1) + x*y
特點方程組：
    -2*lam*x + y
    -2*lam*y + x
    x**2 + y**2 - 1

解：
    {lam: -1/2, x: -sqrt(2)/2, y: sqrt(2)/2}    f = -1/2
    {lam: -1/2, x: sqrt(2)/2, y: -sqrt(2)/2}    f = -1/2
    {lam: 1/2, x: -sqrt(2)/2, y: -sqrt(2)/2}    f = 1/2
    {lam: 1/2, x: sqrt(2)/2, y: sqrt(2)/2}    f = 1/2

=== 案例 2：max x + y  s.t.  x^2 + y^2 = 1 ===
    {lam: -sqrt(2)/2, x: -sqrt(2)/2, y: -sqrt(2)/2}   f = -sqrt(2)
    {lam: sqrt(2)/2, x: sqrt(2)/2, y: sqrt(2)/2}   f = sqrt(2)

=== 案例 3：min x^2 + y^2  s.t.  x + y >= 2 ===
    {lam: 2, x: 1, y: 1}  f = 2

import sympy as sp

x = sp.Symbol('x')
u = sp.Function('u')

print("=== 實驗一：Derivative 物件 ===")
print("sp.Derivative(x**3, x)      =", sp.Derivative(x**3, x))
print("  .doit()                   =", sp.Derivative(x**3, x).doit())
print("sp.diff(u(x), x)            =", sp.diff(u(x), x))
print("sp.diff(u(x), x).doit()     =", sp.diff(u(x), x).doit())
print("未指定變數的 Derivative      =", sp.Derivative(x**3 + x))
print("  .doit()                   =", sp.Derivative(x**3 + x).doit())
print("Derivative 可被當作表達式運算 :", sp.Derivative(x**3, x) + 1)

print()
print("=== 實驗二：微分算子與代入的次序 ===")
expr = x**3 * sp.exp(x)
print("expr            =", expr)
print("sp.diff(expr, x)=", sp.diff(expr, x))
print("diff 後 subs    =", sp.diff(expr, x).subs(x, 2))
print("subs 後 diff    =", sp.diff(expr.subs(x, 2), x))
print("（subs 後 diff 得到 0：因為已經是常數了 —— 這說明 diff 作用在「表達式」上）")

print()
print("=== 實驗三：未求值函數的微分規則 ===")
v = sp.Function('v')
print("d/dx v(x)        =", sp.diff(v(x), x))
print("d/dx v(x)*x      =", sp.diff(v(x) * x, x))
print("d/dx v(x^2)      =", sp.diff(v(x**2), x))
print("d/dx exp(v(x))   =", sp.diff(sp.exp(v(x)), x))
print("d/dx v(sin(x))   =", sp.diff(v(sp.sin(x)), x))

=== 實驗一：Derivative 物件 ===
sp.Derivative(x**3, x)      = Derivative(x**3, x)
  .doit()                   = 3*x**2
sp.diff(u(x), x)            = Derivative(u(x), x)
sp.diff(u(x), x).doit()     = Derivative(u(x), x)
未指定變數的 Derivative      = Derivative(x**3 + x, x)
  .doit()                   = 3*x**2 + 1
Derivative 可被當作表達式運算 : Derivative(x**3, x) + 1

=== 實驗二：微分算子與代入的次序 ===
expr            = x**3*exp(x)
sp.diff(expr, x)= x**3*exp(x) + 3*x**2*exp(x)
diff 後 subs    = 20*exp(2)
subs 後 diff    = 0
（subs 後 diff 得到 0：因為已經是常數了 —— 這說明 diff 作用在「表達式」上）

=== 實驗三：未求值函數的微分規則 ===
d/dx v(x)        = Derivative(v(x), x)
d/dx v(x)*x      = x*Derivative(v(x), x) + v(x)
d/dx v(x^2)      = 2*x*Subs(Derivative(v(_xi_1), _xi_1), _xi_1, x**2)
d/dx exp(v(x))   = exp(v(x))*Derivative(v(x), x)
d/dx v(sin(x))   = cos(x)*Subs(Derivative(v(_xi_1), _xi_1), _xi_1, sin(x))

import numpy as np
import sympy as sp

# ============================================================
# 前向模式 AD：對數數（dual number）
# ============================================================
class Dual:
    """值 v + 一階微分 d * dx。前向模式自動微分的最小實作。"""

    def __init__(self, v, d):
        self.v, self.d = v, d

    @staticmethod
    def _pair(o):
        return (o.v, o.d) if isinstance(o, Dual) else (o, 0.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 __neg__(self):
        return Dual(-self.v, -self.d)

    def __pow__(self, n):
        if isinstance(n, Dual):
            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**n, n * self.v**(n - 1) * self.d)

    def exp(self):
        e = sp.exp(self.v)
        return Dual(e, e * self.d)

    def log(self):
        return Dual(sp.log(self.v), self.d / self.v)

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


# ============================================================
# 反向模式 AD：計算圖 tape + chain rule 反向累積
# ============================================================
class Tape:
    """記錄運算節點，反向時用 chain rule 累積對各輸入的梯度。"""

    def __init__(self):
        self.nodes = []
        self.variables = {}

    def var(self, name, value):
        self.nodes.append((value, (), ()))
        self.variables[name] = len(self.nodes) - 1
        return len(self.nodes) - 1

    def _binop(self, a, b, kind):
        va = self.nodes[a][0]
        vb = self.nodes[b][0]
        if kind == 'add':
            v, g = va + vb, (1.0, 1.0)
        elif kind == 'mul':
            v, g = va * vb, (vb, va)
        else:
            v, g = va / vb, (1.0 / vb, -va / vb**2)
        self.nodes.append((v, (a, b), g))
        return len(self.nodes) - 1

    def add(self, a, b): return self._binop(a, b, 'add')
    def mul(self, a, b): return self._binop(a, b, 'mul')
    def div(self, a, b): return self._binop(a, b, 'div')

    def neg(self, a):
        self.nodes.append((-self.nodes[a][0], (a,), (-1.0,)))
        return len(self.nodes) - 1

    def pow(self, a, n):
        v = self.nodes[a][0]
        self.nodes.append((v**n, (a,), (n * v**(n - 1),)))
        return len(self.nodes) - 1

    def exp(self, a):
        v = self.nodes[a][0]
        e = float(np.exp(v))
        self.nodes.append((e, (a,), (e,)))
        return len(self.nodes) - 1

    def log(self, a):
        v = self.nodes[a][0]
        self.nodes.append((float(np.log(v)), (a,), (1.0 / v,)))
        return len(self.nodes) - 1

    def sin(self, a):
        v = self.nodes[a][0]
        self.nodes.append((float(np.sin(v)), (a,), (float(np.cos(v)),)))
        return len(self.nodes) - 1

    def backward(self, out):
        grads = [0.0] * len(self.nodes)
        grads[out] = 1.0
        for i in range(len(self.nodes) - 1, -1, -1):
            _, parents, local = self.nodes[i]
            g = grads[i]
            for p, lg in zip(parents, local):
                grads[p] += g * lg
        return {name: grads[idx] for name, idx in self.variables.items()}


# ============================================================
# 三方對照：SymPy 符號 / 前向 AD / 反向 AD / 數值梯度
# ============================================================
print("=== 案例 1：f(x) = x^3*exp(-x) + sin(x) ===")
x = sp.Symbol('x')
f_expr = x**3 * sp.exp(-x) + sp.sin(x)
exact = sp.diff(f_expr, x)
x0 = 1.3
print("符號 f'(x)        =", exact)
print("符號 f'(1.3)      =", float(exact.subs(x, x0)))

t = Dual(x0, 1.0)
r = t**3 * (-t).exp() + t.sin()
print("前向 AD 值/導數   =", float(r.v), "/", float(r.d))

tp = Tape()
a = tp.var('x', x0)
d = tp.mul(tp.pow(a, 3), tp.exp(tp.neg(a)))
out = tp.add(d, tp.sin(a))
gr = tp.backward(out)
print("反向 AD 值/導數   =", tp.nodes[out][0], "/", gr['x'])

h = 1e-6
fnum = lambda z: z**3 * np.exp(-z) + np.sin(z)
print("數值梯度（差分）  =", (fnum(x0+h) - fnum(x0-h)) / (2*h))

print()
print("=== 案例 2：f(x,y) = x^2*y + exp(y) - x*y^2 ===")
xa, ya = sp.symbols('x y', real=True)
f2 = xa**2 * ya + sp.exp(ya) - xa * ya**2
gx, gy = sp.diff(f2, xa), sp.diff(f2, ya)
p0 = (0.7, 1.4)
print("符號梯度 =", (gx, gy))
print("符號梯度值 =", (float(gx.subs({xa: p0[0], ya: p0[1]})),
                     float(gy.subs({xa: p0[0], ya: p0[1]}))))

tp2 = Tape()
X = tp2.var('x', p0[0])
Y = tp2.var('y', p0[1])
t1 = tp2.mul(tp2.pow(X, 2), Y)
t3 = tp2.mul(X, tp2.pow(Y, 2))
out2 = tp2.add(tp2.add(t1, tp2.exp(Y)), tp2.neg(t3))
g2 = tp2.backward(out2)
print("反向 AD 梯度 =", (g2['x'], g2['y']))

fx = Dual(p0[0], 1.0); fy0 = Dual(p0[1], 0.0)
vy = fx**2 * fy0 + fy0.exp() - fx * fy0**2
gx2 = Dual(p0[0], 0.0); gy2 = Dual(p0[1], 1.0)
vy2 = gx2**2 * gy2 + gy2.exp() - gx2 * gy2**2
print("前向 AD 梯度 =", (float(vy.d), float(vy2.d)))

hg = 1e-6
f2num = lambda a, b: a**2 * b + np.exp(b) - a * b**2
ndx = (f2num(p0[0]+hg, p0[1]) - f2num(p0[0]-hg, p0[1])) / (2*hg)
ndy = (f2num(p0[0], p0[1]+hg) - f2num(p0[0], p0[1]-hg)) / (2*hg)
print("數值梯度      =", (ndx, ndy))

=== 案例 1：f(x) = x^3*exp(-x) + sin(x) ===
符號 f'(x)        = -x**3*exp(-x) + 3*x**2*exp(-x) + cos(x)
符號 f'(1.3)      = 1.0504826700113055
前向 AD 值/導數   = 1.5623105347129187 / 1.0504826700113055
反向 AD 值/導數   = 1.5623105347129187 / 1.0504826700113055
數值梯度（差分）  = 1.0504826698509717

=== 案例 2：f(x,y) = x^2*y + exp(y) - x*y^2 ===
符號梯度 = (2*x*y - y**2, x**2 - 2*x*y + exp(y))
符號梯度值 = (0.0, 2.585199966844675)
反向 AD 梯度 = (0.0, 2.5851999668446743)
前向 AD 梯度 = (0.0, 2.585199966844675)
數值梯度      = (np.float64(0.0), np.float64(2.585199966542717))

import numpy as np
import sympy as sp

# ============================================================
# 前向 vs 反向：多變數梯度的計算量（以「基本運算次數」計）
# ============================================================
class Tape:
    """最小 tape：支援 +, *, /, neg, pow, exp, log, sin。"""

    def __init__(self):
        self.nodes = []
        self.variables = {}

    def var(self, name, value):
        self.nodes.append((value, (), ()))
        self.variables[name] = len(self.nodes) - 1
        return len(self.nodes) - 1

    def _bin(self, a, b, kind):
        va, vb = self.nodes[a][0], self.nodes[b][0]
        if kind == 'add':
            v, g = va + vb, (1.0, 1.0)
        elif kind == 'mul':
            v, g = va * vb, (vb, va)
        else:
            v, g = va / vb, (1.0 / vb, -va / vb**2)
        self.nodes.append((v, (a, b), g))
        return len(self.nodes) - 1

    def add(self, a, b): return self._bin(a, b, 'add')
    def mul(self, a, b): return self._bin(a, b, 'mul')
    def div(self, a, b): return self._bin(a, b, 'div')

    def neg(self, a):
        self.nodes.append((-self.nodes[a][0], (a,), (-1.0,)))
        return len(self.nodes) - 1

    def pow(self, a, n):
        v = self.nodes[a][0]
        self.nodes.append((v**n, (a,), (n * v**(n - 1),)))
        return len(self.nodes) - 1

    def exp(self, a):
        v = self.nodes[a][0]
        e = float(np.exp(v))
        self.nodes.append((e, (a,), (e,)))
        return len(self.nodes) - 1

    def sin(self, a):
        v = self.nodes[a][0]
        self.nodes.append((float(np.sin(v)), (a,), (float(np.cos(v)),)))
        return len(self.nodes) - 1

    def backward(self, out):
        grads = [0.0]*len(self.nodes)
        grads[out] = 1.0
        for i in range(len(self.nodes)-1, -1, -1):
            _, parents, local = self.nodes[i]
            g = grads[i]
            for p, lg in zip(parents, local):
                grads[p] += g*lg
        return {nm: grads[i] for nm, i in self.variables.items()}

    def count(self):
        """回傳 tape 上的基本運算節點數。"""
        return len(self.nodes)


def forward_cost(n):
    """前向模式：n 個方向，每個方向都跑一條完整的 tape。"""
    total = 0
    for j in range(n):
        tp = Tape()
        outs = []
        for i in range(n):
            xi = tp.var(f"x{i}", 0.0)
            outs.append(tp.mul(tp.sin(xi), tp.exp(xi)))
        acc = outs[0]
        for o in outs[1:]:
            acc = tp.add(acc, o)
        tp.backward(acc)
        total += tp.count()
    return total


def reverse_cost(n):
    """反向模式：1 條 tape，1 次前向 + 1 次反向。"""
    tp = Tape()
    outs = []
    for i in range(n):
        xi = tp.var(f"x{i}", 0.0)
        outs.append(tp.mul(tp.sin(xi), tp.exp(xi)))
    acc = outs[0]
    for o in outs[1:]:
        acc = tp.add(acc, o)
    tp.backward(acc)
    return tp.count()


print("f(x) = sum_i sin(x_i)*exp(x_i)，求全部 n 個偏導")
print(f"{'n':>5} {'前向模式節點數':>16} {'反向模式節點數':>16} {'比值':>8}")
for n in [2, 4, 8, 16, 32]:
    fc = forward_cost(n)
    rc = reverse_cost(n)
    print(f"{n:>5} {fc:>16} {rc:>16} {fc/rc:>7.1f}x")

f(x) = sum_i sin(x_i)*exp(x_i)，求全部 n 個偏導
    n          前向模式節點數          反向模式節點數       比值
    2               18                9     2.0x
    4               76               19     4.0x
    8              312               39     8.0x
   16             1264               79    16.0x
   32             5088              159    32.0x

import numpy as np
import sympy as sp

# ============================================================
# (a) 符號端：推導目標函數與其梯度
# ============================================================
t = sp.symbols('t', real=True)
F = sp.exp(-t**2) * sp.sin(3*t)
F1 = sp.diff(F, t)
F2 = sp.diff(F, t, 2)

print("F(t)  =", F)
print("F'(t) =", F1)
print("F''(t) =", sp.simplify(F2))

# 符號端：求臨界點（如實回報成功或失敗）
try:
    crit = sp.solve(sp.Eq(F1, 0), t)
    print("F' = 0 的解（符號）= ", crit)
    print("解的個數 =", len(crit))
except Exception as e:
    crit = []
    print("F' = 0 的解：符號求解失敗 ->", type(e).__name__)
    print("   （超越方程的符號求解本來就不完整 —— 這是需要數值方法的原因）")

# ============================================================
# (b) AD 端：用 lambdify 把符號梯度編譯成數值函數
# ============================================================
F_num = sp.lambdify(t, F, "numpy")
F1_num = sp.lambdify(t, F1, "numpy")
F2_num = sp.lambdify(t, F2, "numpy")

print()
print("lambdify 後的數值函數：")
print("   F(0.3)  =", F_num(0.3))
print("   F'(0.3) =", F1_num(0.3))
print("   F''(0.3)=", F2_num(0.3))

t0 = 0.30
for k in range(6):
    g = float(F1_num(t0))
    h = float(F2_num(t0))
    t0 = t0 - g/h
    print(f"   Newton step {k}: t = {t0!r}, |F'(t)| = {abs(float(F1_num(t0))):.3e}")

# ============================================================
# (c) 梯度下降：比較「符號梯度」與「數值差分梯度」
# ============================================================
def gradient_descent(grad_fn, x0, lr, steps):
    x = x0
    for _ in range(steps):
        x = x - lr * grad_fn(x)
    return x


def numeric_grad(f, x, h=1e-6):
    return (f(x+h) - f(x-h)) / (2*h)


x_sym = gradient_descent(F1_num, 0.30, 0.1, 30)
x_num = gradient_descent(lambda z: numeric_grad(F_num, z), 0.30, 0.1, 30)
print()
print("用符號梯度做梯度下降 30 步 ->", x_sym)
print("用數值差分做梯度下降 30 步 ->", x_num)
print("兩者之差 =", abs(x_sym - x_num))
print("符號端找到的臨界點 =", crit)

F(t)  = exp(-t**2)*sin(3*t)
F'(t) = -2*t*exp(-t**2)*sin(3*t) + 3*exp(-t**2)*cos(3*t)
F''(t) = (4*t**2*sin(3*t) - 12*t*cos(3*t) - 11*sin(3*t))*exp(-t**2)
F' = 0 的解：符號求解失敗 -> NotImplementedError
   （超越方程的符號求解本來就不完整 —— 這是需要數值方法的原因）

lambdify 後的數值函數：
   F(0.3)  = 0.715906890970694
   F'(0.3) = 1.2747820706516402
   F''(0.3)= -9.662440766209054
   Newton step 0: t = 0.4319316828424694, |F'(t)| = 1.391e-02
   Newton step 1: t = 0.43044541447404705, |F'(t)| = 1.496e-05
   Newton step 2: t = 0.43044701034627586, |F'(t)| = 1.706e-11
   Newton step 3: t = 0.4304470103480949, |F'(t)| = 3.331e-16
   Newton step 4: t = 0.43044701034809485, |F'(t)| = 3.331e-16
   Newton step 5: t = 0.4304470103480949, |F'(t)| = 3.331e-16

用符號梯度做梯度下降 30 步 -> -0.4304470103480949
用數值差分做梯度下降 30 步 -> -0.4304470103525196
兩者之差 = 4.424682842341099e-12
符號端找到的臨界點 = []