3. 自由曲面法向量与曲率计算

任务说明

给定自由曲面方程 z = f(x,y) = 0.3x² + 0.2y² + 0.1xy，编写 Python 程序计算任意点的法向量、主曲率（k₁, k₂）和平均曲率 H，并绘制曲率分布热力图。

数学基础

令 p = ∂z/∂x, q = ∂z/∂y, r = ∂²z/∂x², s = ∂²z/∂x∂y, t = ∂²z/∂y²。

• 法向量（未归一化）：N = (-p, -q, 1)
• 第一基本形式：E=1+p², F=pq, G=1+q²
• 第二基本形式：L=r/√(1+p²+q²), M=s/√(1+p²+q²), N₂=t/√(1+p²+q²)
• 高斯曲率 K=(LN₂-M²)/(EG-F²)，平均曲率 H=(EN₂-2FM+GL)/(2(EG-F²))
• 主曲率 k₁,₂ = H ± √(H²-K)

完整代码

import numpy as np
import matplotlib.pyplot as plt

# 曲面方程 z = 0.3x^2 + 0.2y^2 + 0.1xy
def f(x, y):
    return 0.3*x**2 + 0.2*y**2 + 0.1*x*y

def fx(x, y):   return 0.6*x + 0.1*y
def fy(x, y):   return 0.4*y + 0.1*x
def fxx(x, y):  return np.full_like(x, 0.6)
def fyy(x, y):  return np.full_like(y, 0.4)
def fxy(x, y):  return np.full_like(x, 0.1)

def surface_geometry(x, y):
    p, q = fx(x, y), fy(x, y)
    r, s, t = fxx(x, y), fxy(x, y), fyy(x, y)
    denom = np.sqrt(1 + p**2 + q**2)
    E, F, G = 1 + p**2, p*q, 1 + q**2
    L, M, N2 = r/denom, s/denom, t/denom
    K = (L*N2 - M**2) / (E*G - F**2)
    H = (E*N2 - 2*F*M + G*L) / (2*(E*G - F**2))
    disc = np.sqrt(np.maximum(H**2 - K, 0))
    k1, k2 = H + disc, H - disc
    nvec = np.stack([-p, -q, np.ones_like(p)], axis=-1)
    nvec /= np.linalg.norm(nvec, axis=-1, keepdims=True)
    return nvec, k1, k2, H, K

# 单点验证
x0, y0 = 1.0, 0.5
n, k1, k2, H, K = surface_geometry(np.array(x0), np.array(y0))
print(f"点({x0},{y0}) 法向量={n}")
print(f"k1={k1:.4f}, k2={k2:.4f}, H={H:.4f}, K={K:.4f}")

# 曲率分布可视化
xs = np.linspace(-2, 2, 200)
ys = np.linspace(-2, 2, 200)
X, Y = np.meshgrid(xs, ys)
_, K1, K2, Hm, Km = surface_geometry(X, Y)

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for ax, data, title in zip(axes,
        [K1, K2, Hm],
        ['主曲率 k1', '主曲率 k2', '平均曲率 H']):
    im = ax.contourf(X, Y, data, levels=30, cmap='RdBu_r')
    ax.set_title(title); ax.set_xlabel('x'); ax.set_ylabel('y')
    plt.colorbar(im, ax=ax)
plt.tight_layout()
plt.savefig('curvature_distribution.png', dpi=120)
plt.show()

运行与验证

pip install numpy matplotlib
python curvature.py

验证要点：

• 控制台输出点 (1.0, 0.5) 的法向量、k₁、k₂、H、K 数值。
• 弹出三张子图，分别显示 k₁、k₂、H 在 [-2,2]² 上的分布，并保存为 curvature_distribution.png。
• 通过条件：法向量模长≈1；对抛物面型曲面 H 与 K 应随 |x|、|y| 增大而减小（曲率在原点最大），热力图中心颜色与边缘明显不同。

关键点

代码用解析导数直接计算，避免数值差分误差；np.maximum(H²-K, 0) 防止浮点误差导致根号内出现微小负值。若换成其他曲面方程，只需替换 f 及一阶、二阶偏导函数即可复用整套流程。