一颗心的体积

发表于 分类于 阅读次数：

\[ \left(x^2 + \frac{9}{4}y^2 + z^2 - 1\right)^3 - x^2z^3 - \frac{9}{80}y^2z^3 = 0 \]

十八岁那年，我第一次看到这个曲面方程，被深深地吸引了：简洁的形式、优美的图像。

Taubin Heart Surface

之后的微积分课上，我问老师：怎样计算这个曲面包围的体积呢？ 他说回去思考一下，下次见面告诉我。然后就没有然后了~

后来才知道，这个方程是IBM研究员Gabriel Taubin在1993年构造的。 最初的目的是测试光栅化渲染算法——和浪漫故事无关。

一直记得这颗爱心的形状。 那个搁置了很多年的问题，如今想自己试着回答。 先从证明曲面闭合开始，再用数值方法估算体积。

这颗心完整吗？

数学上要证明该方程描述的是一个闭合曲面，需同时满足闭集、有界与正则性三个条件：

  1. 闭集：连续多项式 \(F(x,y,z)=0\) 的解集（零点原像），拓扑上自然是闭集。
  2. 有界：\(r \to \infty\) 时，\(r^6\) 正项主导 \(r^5\) 余项，在某个有限半径外恒有 \(F > 0\)（方程无解）。
  3. 正则性：通过开方变形消除 \(z=0\) 处的伪奇点，求偏导后利用反证法，可证梯度处处非零。

根据 Heine-Borel 定理，该集合紧致（有界闭集）；由隐函数定理，它是无边界二维流形。两者结合，即为严格的闭合曲面。

心到底有多大？

计算体积需先确定边界。该方程关于 \(z\) 是六次多项式。 Abel-Ruffini 定理指出，五次及以上代数方程不存在通用根式解。 因此边界 \(z(x,y)\) 无法用初等函数表达，体积无初等闭式解。

但是我们可以用数值计算获得近似值。 在 Mathematica 中有两种方法可以交叉验证：

  1. 调用 NIntegrate 对空间区域做数值积分
  2. 生成离散网格，转化为多面体直接求取体积

如下图所示，两种方法的结果在小数点后三位吻合。\(\boxed{V({\color{red}♥}) \approx 3.345}\)

体积计算：积分法与网格法

心里投千亿点

此外，也可以用简单的蒙特卡洛积分来估计心形曲面所围体积。 Julia 实现如下，1000亿采样点得到结果3.3451，4核i5耗时约45s。 代码使用SplitMix64生成随机数，避开了多线程下RNG的锁竞争。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
using Base.Threads

@inline @fastmath function is_inside(x, y, z)
    ((x^2 + 2.25y^2 + z^2 - 1.0)^3 - (x^2 + 0.1125y^2) * z^3) < 0
end

@inline function next_rand(s::UInt64)
    s += 0x9e3779b97f4a7c15
    z = xor(s, s >> 30) * 0xbf58476d1ce4e5b9
    z = xor(z, z >> 27) * 0x94d049bb133111eb
    u01 = Float64(xor(z, z >> 31) >> 11) * 2.0^-53
    return u01, s
end

function mc_vol(n, nthr = Threads.nthreads())
    chunk, remainder = divrem(n, nthr)
    hits = zeros(Int64, nthr)
    @fastmath Threads.@threads for tid in 1:nthr
        samples = chunk + (tid <= remainder ? 1 : 0)
        state, count = UInt64(12345 + tid * 31337), 0
        for _ in 1:samples
            ux, state = next_rand(state)
            uy, state = next_rand(state)
            uz, state = next_rand(state)
            count += is_inside(3ux - 1.5, 3uy - 1.5, 3uz - 1.5)
        end
        @inbounds hits[tid] = count
    end
    sum(hits)
end

mc_vol(10_000)  # warmup
let N = 100_000_000_000
    elapsed = @elapsed hits = mc_vol(N)
    volume = (hits / N) * 27.0
    println("Volume: $(round(volume, digits=4))")
    println("Time:   $(round(elapsed, digits=3))s")
end

最后的碎碎念

那个画心形的少年，和今晚的中年大叔，是同一个人。 当年没等到答案的问题，如今求出了近似值。

可有些东西，就像心形曲面的精确体积，存在、可逼近、无法到达。 今日5.21，农历小满。人生小满胜万全。

参考资料

G. Taubin. Rasterizing Algebraic Curves and Surfaces. IEEE Computer Graphics and Applications, 1994, 14(2): 14-23.