9.3 用 SymPy 計算曲線曲率與曲面度量#

import sympy as sp

def clean(expr):
    """simplify 對三角函數有時不夠強，補一層 trigsimp。"""
    return sp.simplify(sp.trigsimp(sp.simplify(expr)))

def tidy(expr, t):
    """把三角式整理成以 cos(t) 為主的形式（給輸出報表用）。"""
    e = sp.simplify(sp.expand_trig(expr))
    e = e.subs(sp.sin(t)**2, 1 - sp.cos(t)**2)
    return sp.simplify(sp.expand_trig(e))

def frenet_frame(r, t):
    """參數化曲線 r(t) (SymPy Matrix) -> (T, N, B, kappa, tau)。"""
    r1 = sp.Matrix([clean(x) for x in r.diff(t)])
    spd = clean(r1.norm())                        # 速率 |r'|
    T = sp.Matrix([clean(x) for x in r1 / spd])    # 單位切向量
    T1 = sp.Matrix([clean(x) for x in T.diff(t)])
    T1n = clean(T1.norm())
    kappa = clean(T1n / spd)                      # kappa = |dT/ds|
    N = sp.Matrix([clean(x) for x in T1 / T1n])
    B = sp.Matrix([clean(x) for x in T.cross(N)])
    tau = clean(-sp.Matrix([clean(x) for x in B.diff(t)]).dot(N) / spd)
    return T, N, B, kappa, tau

def plane_kappa(xf, yf, t):
    """平面曲線 (x(t), y(t)) 的曲率 |x'y''-y'x''| / (x'^2+y'^2)^(3/2)。"""
    r = sp.Matrix([xf, yf, sp.Integer(0)])
    r1 = r.diff(t); cr = r1.cross(r.diff(t, 2))
    return clean(cr.norm() / r1.norm()**3)

def surface_forms(r, u, v, normal_manual=None):
    """參數曲面 r(u,v) -> (E,F,G), (L,M,N2), det_g, K, H, n, g, b。
    記號取捨：第二基本形式的第三個分量記作 N2（避免與法向量 n 混淆）。"""
    ru = r.diff(u); rv = r.diff(v)
    E = clean(ru.dot(ru)); F = clean(ru.dot(rv)); G = clean(rv.dot(rv))
    if normal_manual is None:
        cru = clean(ru.cross(rv))
        n = sp.Matrix([clean(x) for x in cru / cru.norm()])
    else:
        n = normal_manual
    L = clean(r.diff(u, 2).dot(n))        # r_uu . n
    M = clean(r.diff(u, 1, v, 1).dot(n))  # r_uv . n
    N2 = clean(r.diff(v, 2).dot(n))       # r_vv . n
    det_g = clean(E * G - F**2)
    return dict(E=E, F=F, G=G, L=L, M=M, N2=N2, det_g=det_g,
                K=clean((L*N2 - M**2) / det_g),
                H=clean((E*N2 - 2*F*M + G*L) / (2 * det_g)),
                n=n, g=sp.Matrix([[E, F], [F, G]]), b=sp.Matrix([[L, M], [M, N2]]))

def christoffel(E, F, G, u, v):
    """第一基本形式 -> Christoffel 符號 Gamma[m][A][B]（0<->u, 1<->v）。"""
    g = sp.Matrix([[E, F], [F, G]]); ginv = g.inv(); coords = [u, v]
    gd = [[[sp.diff(g[i, j], coords[l]) for l in range(2)]
            for j in range(2)] for i in range(2)]
    Gamma = [[[None]*2 for _ in range(2)] for _ in range(2)]
    for m in range(2):
        for A in range(2):
            for B in range(2):
                Gamma[m][A][B] = sp.simplify(sum(
                    ginv[m, c]*(gd[B][c][A] + gd[A][c][B] - gd[A][B][c])/2
                    for c in range(2)))
    return Gamma

def intrinsic_K(E, F, G, u, v):
    """內蘊高斯曲率 K = R_{1212} / (EG - F^2)：只用 E, F, G。"""
    Gamma = christoffel(E, F, G, u, v); coords = [u, v]
    def ri(r, sidx, mu, nu):      # R^r_{s mu nu}，符號約定見 9.1
        s = (sp.diff(Gamma[r][nu][sidx], coords[mu])
             - sp.diff(Gamma[r][mu][sidx], coords[nu]))
        for lam in range(2):
            s += (Gamma[r][mu][lam]*Gamma[lam][nu][sidx]
                  - Gamma[r][nu][lam]*Gamma[lam][mu][sidx])
        return sp.simplify(s)
    R1212 = sp.simplify(E*ri(0, 1, 0, 1) + F*ri(1, 1, 0, 1))
    return sp.simplify(R1212 / (E*G - F**2))

def geodesic_residuals(E, F, G, u, v):
    """測地線方程左邊（殘差）：(eq_u, eq_v)，其中 u = us(s), v = vs(s)。"""
    Gamma = christoffel(E, F, G, u, v)
    s = sp.Symbol('s', real=True)
    us = sp.Function('us')(s); vs = sp.Function('vs')(s)
    eq_u = sp.simplify(sp.diff(us, s, 2)
                        + Gamma[0][0][1]*sp.diff(us, s)*sp.diff(vs, s)
                        + Gamma[0][1][1]*sp.diff(vs, s)**2)
    eq_v = sp.simplify(sp.diff(vs, s, 2)
                        + Gamma[1][0][0]*sp.diff(us, s)**2
                        + 2*Gamma[1][0][1]*sp.diff(us, s)*sp.diff(vs, s))
    return eq_u, eq_v

def check_geodesic(eq_u, eq_v, u, uval, vval):
    """把 (u(s), v(s)) 的解代進測地線方程，檢查殘差是否為 0。"""
    s = sp.Symbol('s', real=True)
    subs = {sp.Function('us')(s): uval, sp.Function('vs')(s): vval}
    return (sp.simplify(eq_u.subs(subs).subs(u, uval)),
            sp.simplify(eq_v.subs(subs).subs(u, uval)))

import sympy as sp
from geo_lib import *

t = sp.symbols('t', real=True)
a, b = sp.symbols('a b', positive=True)

r_helix = sp.Matrix([a*sp.cos(t), a*sp.sin(t), b*t])
T, N, B, kappa, tau = frenet_frame(r_helix, t)
print("T =", T.T); print("N =", N.T); print("B =", B.T)
print("kappa =", kappa, " tau =", tau)
print("正交性 T.N, T.B, N.B =", clean(T.dot(N)), clean(T.dot(B)), clean(N.dot(B)))
print("單位性 |T|,|N|,|B| =", clean(T.norm()), clean(N.norm()), clean(B.norm()))

# 驗證 Frenet–Serret 方程：dT/ds = kappa N, dN/ds = -kappa T + tau B, dB/ds = -tau N
r1 = r_helix.diff(t); spd = clean(r1.norm())
def dds(V): return sp.Matrix([clean(x) for x in V.diff(t)]) / spd
for label, residual in [
        ("dT/ds - kappa*N",        dds(T) - kappa*N),
        ("dN/ds - (-kappa T + tau B)", dds(N) - (-kappa*T + tau*B)),
        ("dB/ds + tau*N",          dds(B) + tau*N)]:
    print(label, "=", sp.Matrix([clean(sp.simplify(residual[i])) for i in range(3)]))

# 用 9.2 的行列式公式交叉驗證 kappa 與 tau
cr = r_helix.diff(t).cross(r_helix.diff(t, 2))
print("kappa (cross) =", clean(cr.norm()/spd**3), " 差 =", clean(kappa - cr.norm()/spd**3))
print("tau   (det)   =", clean(cr.dot(r_helix.diff(t, 3))/cr.dot(cr)),
      " 差 =", clean(tau - cr.dot(r_helix.diff(t, 3))/cr.dot(cr)))
print("arc length 0->s =", sp.simplify(sp.integrate(spd, (t, 0, sp.Symbol('s', real=True)))))
# ==>
# T = Matrix([[-a*sin(t)/sqrt(a**2 + b**2), a*cos(t)/sqrt(a**2 + b**2), b/sqrt(a**2 + b**2)]])
# N = Matrix([[-cos(t), -sin(t), 0]])
# B = Matrix([[b*sin(t)/sqrt(a**2 + b**2), -b*cos(t)/sqrt(a**2 + b**2), a/sqrt(a**2 + b**2)]])
# kappa = a/(a**2 + b**2)  tau = b/(a**2 + b**2)
# 正交性 T.N, T.B, N.B = 0 0 0
# 單位性 |T|,|N|,|B| = 1 1 1
# dT/ds - kappa*N = Matrix([[0], [0], [0]])
# dN/ds - (-kappa T + tau B) = Matrix([[0], [0], [0]])
# dB/ds + tau*N = Matrix([[0], [0], [0]])
# kappa (cross) = a/(a**2 + b**2)  差 = 0
# tau   (det)   = b/(a**2 + b**2)  差 = 0
# arc length 0->s = s*sqrt(a**2 + b**2)

r_cat = sp.Matrix([t, a*sp.cosh(t/a), sp.Integer(0)])
T, N, B, kappa, tau = frenet_frame(r_cat, t)
print("懸鏈線 T =", T.T); print("懸鏈線 N =", N.T); print("懸鏈線 B =", B.T)
print("懸鏈線 kappa =", kappa, " tau =", tau, " kappa at t=0 =", kappa.subs(t, 0))

r_viv = sp.Matrix([a*sp.cos(t)**2, a*sp.sin(t)*sp.cos(t), a*sp.sin(t)])
print("Viviani 在球面上 |r|^2 =", sp.simplify(r_viv.dot(r_viv)))
print("Viviani 在圓柱上 x^2+y^2 - a*x =",
      sp.simplify(r_viv[0]**2 + r_viv[1]**2 - a*r_viv[0]))
T, N, B, kappa, tau = frenet_frame(r_viv, t)
print("Viviani kappa =", tidy(kappa, t)); print("Viviani tau =", tidy(tau, t))
# ==>
# 懸鏈線 T = Matrix([[1/cosh(t/a), tanh(t/a), 0]])
# 懸鏈線 N = Matrix([[-tanh(t/a), 1/cosh(t/a), 0]])
# 懸鏈線 B = Matrix([[0, 0, 1]])
# 懸鏈線 kappa = 1/(a*cosh(t/a)**2)  tau = 0  kappa at t=0 = 1/a
# Viviani 在球面上 |r|^2 = a**2
# Viviani 在圓柱上 x^2+y^2 - a*x = 0
# Viviani kappa = sqrt(3*cos(t)**2 + 5)/(a*(cos(t)**2 + 1)**(3/2))
# Viviani tau = 6*cos(t)/(a*(3*cos(t)**2 + 5))

import numpy as np

def frenet_numeric(r, t, t0, aval, h=1e-6):
    """中心差分獨立算出 (kappa, tau)，作為符號結果的對照組。"""
    aS = sp.Symbol('a', positive=True)
    f1 = sp.lambdify((t, aS), sp.diff(r, t), 'numpy')
    def der(x): return np.asarray(f1(x, aval), dtype=float).ravel()
    def nrm(v): return float(np.sqrt(np.sum(v**2)))
    def T(x): v = der(x); return v/nrm(v)
    def N(x):
        dT = (T(x+h)-T(x-h))/(2*h); return dT/nrm(dT)
    def B(x): return np.cross(T(x), N(x))
    dT = (T(t0+h)-T(t0-h))/(2*h); dB = (B(t0+h)-B(t0-h))/(2*h)
    speed = nrm(der(t0))
    return nrm(dT)/speed, -float(np.dot(dB, N(t0)))/speed

kappa_cand = sp.sqrt(3*sp.cos(t)**2 + 5)/(a*(sp.cos(t)**2 + 1)**sp.Rational(3, 2))
tau_cand = 6*sp.cos(t)/(a*(3*sp.cos(t)**2 + 5))
for t0 in [0.0, 0.3, 0.7, 1.0, 1.5]:
    kn, tn = frenet_numeric(r_viv, t, t0, 1.0)
    print(f"t={t0}: kappa num={kn:.8f} sym={float(kappa_cand.subs({t:t0, a:1.0})):.8f}"
          f" | tau num={tn:.8f} sym={float(tau_cand.subs({t:t0, a:1.0})):.8f}")
# ==>
# t=0.0: kappa num=1.00000000 sym=1.00000000 | tau num=0.74999392 sym=0.75000000
# t=0.3: kappa num=1.05161104 sym=1.05161104 | tau num=0.74076010 sym=0.74076201
# t=0.7: kappa num=1.30248864 sym=1.30248864 | tau num=0.67937766 sym=0.67936145
# t=1.0: kappa num=1.65072984 sym=1.65072984 | tau num=0.55172707 sym=0.55172487
# t=1.5: kappa num=2.22271833 sym=2.22271833 | tau num=0.08463017 sym=0.08463056

import sympy as sp
from geo_lib import *
import numpy as np
import matplotlib.pyplot as plt

t = sp.symbols('t', real=True)
a, b = sp.symbols('a b', positive=True)

X = sp.Function('X')(t); Y = sp.Function('Y')(t)
r = sp.Matrix([X, Y, sp.Integer(0)])
r1 = r.diff(t); cross = r1.cross(r.diff(t, 2))
print("r' x r'' =", cross.T)
print("|r'|^2   =", sp.simplify(r1.dot(r1)))
print("kappa    =", sp.simplify(sp.sqrt(cross.dot(cross))/(sp.sqrt(r1.dot(r1)))**3))

print("圓     k =", plane_kappa(a*sp.cos(t), a*sp.sin(t), t))
print("橢圓   k =", plane_kappa(a*sp.cos(t), b*sp.sin(t), t))
print("拋物線 k =", plane_kappa(t, t**2, t))
print("雙曲線 k =", plane_kappa(a*sp.cosh(t), b*sp.sinh(t), t))
print("懸鏈線 k =", plane_kappa(t, a*sp.cosh(t/a), t))

ts = np.linspace(-4, 4, 400)
f_par = sp.lambdify(t, plane_kappa(t, t**2, t), 'numpy')
f_hyp = sp.lambdify((t, a, b), plane_kappa(a*sp.cosh(t), b*sp.sinh(t), t), 'numpy')
fig, ax = plt.subplots(1, 2, figsize=(9, 3.4))
ax[0].plot(ts, f_par(ts)); ax[1].plot(ts, f_hyp(ts, 2.0, 1.0), color='C2')
ax[0].set_title(r'parabola $y=x^2$: $\kappa=2/(4t^2+1)^{3/2}$')
ax[1].set_title(r'hyperbola $a=2,b=1$')
for A in ax: A.set_xlabel('t'); A.set_ylabel(r'$\kappa$')
plt.tight_layout(); plt.savefig('/tmp/geo93_curv.png', dpi=90)
print("kappa parabola at t=0 :", float(f_par(0.0)))
# ==>
# r' x r'' = Matrix([[0, 0, Derivative(X(t), t)*Derivative(Y(t), (t, 2)) - Derivative(X(t), (t, 2))*Derivative(Y(t), t)]])
# |r'|^2   = Derivative(X(t), t)**2 + Derivative(Y(t), t)**2
# kappa    = sqrt((Derivative(X(t), t)*Derivative(Y(t), (t, 2)) - Derivative(X(t), (t, 2))*Derivative(Y(t), t))**2)/(Derivative(X(t), t)**2 + Derivative(Y(t), t)**2)**(3/2)
# 圓     k = 1/a
# 橢圓   k = a*b/(a**2*sin(t)**2 + b**2*cos(t)**2)**(3/2)
# 拋物線 k = 2/(4*t**2 + 1)**(3/2)
# 雙曲線 k = a*b/(a**2*sinh(t)**2 + b**2*cosh(t)**2)**(3/2)
# 懸鏈線 k = 1/(a*cosh(t/a)**2)
# kappa parabola at t=0 : 2.0

import sympy as sp
from geo_lib import *

R = sp.symbols('R', positive=True)
u, v = sp.symbols('u v', real=True)
a = sp.symbols('a', positive=True)
x, y = sp.symbols('x y', real=True)

r_sph = sp.Matrix([R*sp.sin(u)*sp.cos(v), R*sp.sin(u)*sp.sin(v), R*sp.cos(u)])
# 用幾何直覺指定法向量：朝外（= r/R），避開 Abs(sin u) 噪音
n_out = sp.Matrix([sp.sin(u)*sp.cos(v), sp.sin(u)*sp.sin(v), sp.cos(u)])
d = surface_forms(r_sph, u, v, normal_manual=n_out)
print("球面 E =", d['E'], "| F =", d['F'], "| G =", d['G'])
print("球面 L =", d['L'], "| M =", d['M'], "| N2 =", d['N2'])
print("球面 det g =", d['det_g'], "| K =", d['K'], "| H =", d['H'])
print("g 矩陣 =\n", d['g']); print("b 矩陣 =\n", d['b'])

r_cyl = sp.Matrix([a*sp.cos(v), a*sp.sin(v), u])
d = surface_forms(r_cyl, u, v)
print("圓柱 E,F,G =", d['E'], d['F'], d['G'])
print("圓柱 L,M,N2 =", d['L'], d['M'], d['N2'])
print("圓柱 K =", d['K'], "| H =", d['H'], "| det g =", d['det_g'])

r_sad = sp.Matrix([x, y, x**2 - y**2])
d = surface_forms(r_sad, x, y)
print("鞍面 E,F,G =", d['E'], d['F'], d['G'])
print("鞍面 L,M,N2 =", d['L'], d['M'], d['N2'])
print("鞍面 K =", d['K'], "| H =", d['H'])
print("K at (0,0) =", d['K'].subs({x: 0, y: 0}))
print("K 與 -4/(1+4x^2+4y^2)^2 的差 =",
      sp.simplify(d['K'] - (-4/(1 + 4*x**2 + 4*y**2)**2)))
# ==>
# 球面 E = R**2 | F = 0 | G = R**2*sin(u)**2
# 球面 L = -R | M = 0 | N2 = -R*sin(u)**2
# 球面 det g = R**4*sin(u)**2 | K = R**(-2) | H = -1/R
# g 矩陣 =
#  Matrix([[R**2, 0], [0, R**2*sin(u)**2]])
# b 矩陣 =
#  Matrix([[-R, 0], [0, -R*sin(u)**2]])
# 圓柱 E,F,G = 1 0 a**2
# 圓柱 L,M,N2 = 0 0 a
# 圓柱 K = 0 | H = 1/(2*a) | det g = a**2
# 鞍面 E,F,G = 4*x**2 + 1 -4*x*y 4*y**2 + 1
# 鞍面 L,M,N2 = 2/sqrt(4*x**2 + 4*y**2 + 1) 0 -2/sqrt(4*x**2 + 4*y**2 + 1)
# 鞍面 K = -4/(4*x**2 + 4*y**2 + 1)**2 | H = 4*(-x**2 + y**2)/(4*x**2 + 4*y**2 + 1)**(3/2)
# K at (0,0) = -4
# K 與 -4/(1+4x^2+4y^2)^2 的差 = 0

import sympy as sp
from geo_lib import *

R = sp.symbols('R', positive=True)
u, v = sp.symbols('u v', real=True)
r_sph = sp.Matrix([R*sp.sin(u)*sp.cos(v), R*sp.sin(u)*sp.sin(v), R*sp.cos(u)])
n_out = sp.Matrix([sp.sin(u)*sp.cos(v), sp.sin(u)*sp.sin(v), sp.cos(u)])
d = surface_forms(r_sph, u, v, normal_manual=n_out)

print("正確 K =", d['K'])
L_wrong = d['L']/2                       # 誤把 Taylor 展開的 1/2 塞進 L
K_wrong = sp.simplify((L_wrong*d['N2'] - d['M']**2)/d['det_g'])
print("誤算 K =", K_wrong)
print("誤算與正確的比值 =", sp.simplify(K_wrong/d['K']))

dA = sp.sqrt(d['det_g'])
total = sp.simplify(sp.integrate(sp.integrate(d['K']*dA, (u, 0, sp.pi)), (v, 0, 2*sp.pi)))
area = sp.simplify(sp.integrate(sp.integrate(dA, (u, 0, sp.pi)), (v, 0, 2*sp.pi)))
print("總曲率 =", total, " 面積 =", area, " 總曲率/面積 =", sp.simplify(total/area))
# ==>
# 正確 K = R**(-2)
# 誤算 K = 1/(2*R**2)
# 誤算與正確的比值 = 1/2
# 總曲率 = 4*pi 面積 = 4*pi*R**2 總曲率/面積 = R**(-2)

import sympy as sp
from geo_lib import *

u, v = sp.symbols('u v', real=True)
R = sp.symbols('R', positive=True)
a = sp.symbols('a', positive=True)
f = sp.Function('f')
fp = sp.diff(f(u), u); fpp = sp.diff(f(u), u, 2)
# 母線切向量 r_u = (f' cos v, f' sin v, 1)，繞 z 軸旋轉
# 法向量取 (-cos v, -sin v, f') / sqrt(1+f'^2)
n_rev = sp.Matrix([-sp.cos(v)/sp.sqrt(1 + fp**2),
                   -sp.sin(v)/sp.sqrt(1 + fp**2),
                   fp/sp.sqrt(1 + fp**2)])
r_rev = sp.Matrix([f(u)*sp.cos(v), f(u)*sp.sin(v), u])

d = surface_forms(r_rev, u, v, normal_manual=n_rev)
print("E =", d['E'], "| F =", d['F'], "| G =", d['G'])
print("L =", d['L'], "| M =", d['M'], "| N2 =", d['N2'])
print("K =", d['K'])
target = -fpp/(f(u)*(1 + fp**2)**2)
print("目標式 -f''/(f(1+f'^2)^2) =", target)
print("差 =", sp.simplify(d['K'] - target))

# 拿公式去算實際曲面：catenoid 與球面
def rev_normal(ff, u, v):
    fp = sp.diff(ff, u)
    return sp.Matrix([-sp.cos(v)/sp.sqrt(1 + fp**2),
                      -sp.sin(v)/sp.sqrt(1 + fp**2),
                      fp/sp.sqrt(1 + fp**2)])

f_cat = a*sp.cosh(u/a)
r_cat = sp.Matrix([f_cat*sp.cos(v), f_cat*sp.sin(v), u])
d_cat = surface_forms(r_cat, u, v, normal_manual=rev_normal(f_cat, u, v))
print("catenoid K =", sp.simplify(sp.trigsimp(d_cat['K'])))
print("= -1/(a^2 cosh^4(u/a)) ? 差 =",
      sp.simplify(sp.trigsimp(d_cat['K']) + 1/(a**2*sp.cosh(u/a)**4)))

f_s = sp.sqrt(R**2 - u**2)                # 球面（以 z 為高度參數）
r_s = sp.Matrix([f_s*sp.cos(v), f_s*sp.sin(v), u])
d_s = surface_forms(r_s, u, v, normal_manual=rev_normal(f_s, u, v))
print("球面（另一參數化）K =", sp.simplify(sp.radsimp(d_s['K'])))
# ==>
# E = Derivative(f(u), u)**2 + 1 | F = 0 | G = f(u)**2
# L = -Derivative(f(u), (u, 2))/sqrt(Derivative(f(u), u)**2 + 1) | M = 0 | N2 = f(u)/sqrt(Derivative(f(u), u)**2 + 1)
# K = -Derivative(f(u), (u, 2))/((Derivative(f(u), u)**2 + 1)**2*f(u))
# 目標式 -f''/(f(1+f'^2)^2) = -Derivative(f(u), (u, 2))/((Derivative(f(u), u)**2 + 1)**2*f(u))
# 差 = 0
# catenoid K = -1/(a**2*cosh(u/a)**4)
# = -1/(a^2 cosh^4(u/a)) ? 差 = 0
# 球面（另一參數化）K = R**(-2)

import sympy as sp
from geo_lib import *

R = sp.symbols('R', positive=True)
u, v = sp.symbols('u v', real=True)
a = sp.symbols('a', positive=True)
x, y = sp.symbols('x y', real=True)
cases = []

r_sph = sp.Matrix([R*sp.sin(u)*sp.cos(v), R*sp.sin(u)*sp.sin(v), R*sp.cos(u)])
n_out = sp.Matrix([sp.sin(u)*sp.cos(v), sp.sin(u)*sp.sin(v), sp.cos(u)])
cases.append(("球面", surface_forms(r_sph, u, v, normal_manual=n_out), u, v))

r_cyl = sp.Matrix([a*sp.cos(v), a*sp.sin(v), u])
cases.append(("圓柱面", surface_forms(r_cyl, u, v), u, v))

r_sad = sp.Matrix([x, y, x**2 - y**2])
cases.append(("鞍面", surface_forms(r_sad, x, y), x, y))

f = sp.Function('f'); fp = sp.diff(f(u), u)
n_rev = sp.Matrix([-sp.cos(v)/sp.sqrt(1 + fp**2),
                   -sp.sin(v)/sp.sqrt(1 + fp**2),
                   fp/sp.sqrt(1 + fp**2)])
r_rev = sp.Matrix([f(u)*sp.cos(v), f(u)*sp.sin(v), u])
cases.append(("旋轉曲面", surface_forms(r_rev, u, v, normal_manual=n_rev), u, v))

for name, d, cu, cv in cases:
    K_ext = d['K']; K_int = intrinsic_K(d['E'], d['F'], d['G'], cu, cv)
    print(f"--- {name} ---")
    print("  外蘊 K =", K_ext)
    print("  內蘊 K =", K_int)
    print("  差     =", sp.simplify(K_ext - K_int))
# ==>
# --- 球面 ---
#   外蘊 K = R**(-2)
#   內蘊 K = R**(-2)
#   差     = 0
# --- 圓柱面 ---
#   外蘊 K = 0
#   內蘊 K = 0
#   差     = 0
# --- 鞍面 ---
#   外蘊 K = -4/(4*x**2 + 4*y**2 + 1)**2
#   內蘊 K = -4/(16*x**4 + 32*x**2*y**2 + 8*x**2 + 16*y**4 + 8*y**2 + 1)
#   差     = 0
# --- 旋轉曲面 ---
#   外蘊 K = -Derivative(f(u), (u, 2))/((Derivative(f(u), u)**2 + 1)**2*f(u))
#   內蘊 K = -Derivative(f(u), (u, 2))/((Derivative(f(u), u)**2 + 1)**2*f(u))
#   差     = 0

import sympy as sp
from geo_lib import *

R = sp.symbols('R', positive=True)
u, v = sp.symbols('u v', real=True)
a = sp.symbols('a', positive=True)
Gamma = christoffel(R**2, 0, R**2*sp.sin(u)**2, u, v)
names = ['u', 'v']
for m, A, B in [(m, A, B) for m in range(2) for A in range(2) for B in range(2)]:
    if Gamma[m][A][B] != 0:
        print(f"Gamma^{names[m]}_{names[A]}{names[B]} =", Gamma[m][A][B])

eq_u, eq_v = geodesic_residuals(R**2, 0, R**2*sp.sin(u)**2, u, v)
print("eq_u =", eq_u); print("eq_v =", eq_v)

s = sp.Symbol('s', real=True); v_0 = sp.Symbol('v_0')
print("赤道 殘差 =", check_geodesic(eq_u, eq_v, u, sp.pi/2, s/R))   # u = pi/2, v = s/R
print("子午線 殘差 =", check_geodesic(eq_u, eq_v, u, s/R, v_0))        # u = s/R, v = v_0

eq_u, eq_v = geodesic_residuals(1, 0, a**2, u, v)
print("圓柱 eq_u =", eq_u, " eq_v =", eq_v)
print("螺旋線 u=s, v=2s 殘差 =", check_geodesic(eq_u, eq_v, u, s, 2*s))
# ==>
# Gamma^u_vv = -sin(2*u)/2
# Gamma^v_uv = 1/tan(u)
# Gamma^v_vu = 1/tan(u)
# eq_u = -sin(2*u)*Derivative(vs(s), s)**2/2 + Derivative(us(s), (s, 2))
# eq_v = Derivative(vs(s), (s, 2)) + 2*Derivative(us(s), s)*Derivative(vs(s), s)/tan(u)
# 赤道 殘差 = (0, 0)
# 子午線 殘差 = (0, 0)
# 圓柱 eq_u = Derivative(us(s), (s, 2))  eq_v = Derivative(vs(s), (s, 2))
# 螺旋線 u=s, v=2s 殘差 = (0, 0)