PDE，广义特征问题，和神经特征函数法

为了读今年的best paper我也是拼了，天知道我本科只学过ODE，还全是靠背诵考试的（倒

paper的题目是Shape Space Spectra(Yue Chang, Otman Benchekroun, Maurizio M. Chiaramonte, Peter Yichen Chen, and Eitan Grinspun. 2025. Shape Space Spectra. ACM Trans. Graph. 44, 4, Article 121 (August 2025), 16 pages. https://doi.org/10.1145/3731148)，博客会从（我作为cs专业本科没学过的）物理和数学基础开始写，一直写到能看懂这篇文章是怎么做的。

PDE表示物理问题

大部分的物理问题，可以表示为二阶偏微分方程的形式，例如\(F=ma\)，a自然包含了\(u(x, t)\)（设为位移场）对t的二阶导数。比如热传导的问题，\(u_t(x, t) = u_{xx}(x, t)\)，自然包含了\(u(x, t)\)对x的二阶导。等等等等。查gpt会说，物理公式背后一般三种定律，守恒律（conservation laws），力学/几何原理，波动/谐振问题，都和场二阶梯度有关；所以二阶导（拉普拉斯算子）是很常见的公式核心。有一些不包括在内的，比如波动方程，流体力学，这里不讨论。

对于这类问题，在边界条件不同的时候，即使是一样的物理定律，得到的\(u(x, t)\)方程也会不一样，这很好理解数学上ODE也是如此，物理上对不同的形状（形状边界不同）的弹性物体施加同样的力它们的运动是不同的，对同样形状的弹性物体施加不同的力（初始值边界不同）它们的运动也是不同的，尽管背后的物理规律PDE是完全相同的。

这导致这类问题对于计算机来说会很难解，不可能每次都数学去推\(u(x, t)\)是什么。数值的解法，是用离散的时间步，每一步PDE和数值微分公式，可以连立一个线性方程组，每一步解方程组（求矩阵的逆）求解下一步的状态。这个速度无疑是很慢的。

在进入下一部分特征函数解法之前，先记录一下刚度矩阵（K）和质量矩阵（M），我理解这种范式是为了计算机的离散仿真方便做的，在不少论文里面其实会直接跳过物理连续的PDE公式（强形式）直接写基于刚度矩阵和质量矩阵的离散形式，且这个形式本身也很容易理解物理意义。这里以线性弹性动力学（也是SSS这篇文章考虑的范围）为例推一下，实际其他各种问题都可以表示成KM的离散形式。

\[\rho \ddot{u}(x, t)=\nabla \cdot \sigma(u(x, t)) \]

这里\(\rho\)是密度，\(u(x, t)\)是位移场，\(\ddot{u}(x, t)\)是对t的二阶导也就是加速度，所以左边其实就是ma；\(\sigma=C: \varepsilon(u)\)，\(C\)是弹性系数，\(\varepsilon(u)=\frac{1}{2}\left(\nabla u+\nabla u^T\right.)\)被称为应变（strain）是局部拉伸/剪切形变，\(\sigma\)被称为应力（stress）是局部拉伸/剪切形变产生的力。应力是一个3x3矩阵，表示总的应力分解到不同方向的强度，对应应变是各个不同方向的形变。这么做的意义我理解是因为C是不同方向不同的。然后通过散度\(\nabla\)映射回3x1向量也就是xyz三个方向的力，向量力。所有右边其实是F，这个式子也可以理解为F=ma。

对式子两边都乘一个任意函数（被称为测试函数）v，并在物体形状域\(\Omega\)积分，这一步可以理解为整个物体做任意位移，则每个形状位置都做虚功，等式两边仍然成立。

\[\int_{\Omega} v \cdot (\nabla \cdot \sigma(u(x, t)) - \rho \ddot{u}(x, t)) dx=0 \]

拆开之后，对散度项做分部积分，具体来说，对于\(\int_{\Omega} v \cdot (\nabla\cdot \sigma)dx\)，以其中一项为例：

\[\int_\Omega v_i \frac{\partial\sigma_{ij}}{\partial x_j} dx= \int_\Omega \frac{\partial(\sigma_{ij}v_i)}{\partial x_j} dx- \int_\Omega \frac{\partial v_i}{\partial x_j} \sigma_j dx \]

这一步是典型的分部积分，接下来对第一项考虑高斯散度定理，具体来说，体积积分的散度=边界积分。可以理解为比如有一块空气，设定虚拟的边界，对内部气流微分的积分，就等于计算总的流入流出的气流，也就等于对边界的积分。因此可以替换第一项为\(\int_{\partial \Omega} v_i \sigma_{ij} n_j dx\)，这里n是边界法向量。边界项全部暂时不考虑，总体写成：

\[\int_{\Omega} \rho \ddot{u} \cdot v d x+\int_{\Omega} \varepsilon(v): C: \varepsilon(u) d x=\text { (boundary terms). } \]

这个形式被称为弱形式。

接下来做离散化，\(u=NU, v=NV\)，U是控制点被称为节点位移向量，比如说可以理解为mash的顶点，N是把U插值回任意(x,y)的插值函数，被称为形函数矩阵。注意，离散化后，只有N依赖于x，UV不再依赖于x。

则有\(\varepsilon = BU\)，B是N的导数；\(\sigma=D\varepsilon\)，D是材料系数。

那么上面的两项可以写成：

$第一项=\int_{\Omega} \rho \ddot{u} \cdot v d x = \int (NV)^T \rho (N\ddot{U}) dx = V^T (\int N^T\rho N dx) \ddot{U} = V^TM\ddot{U} $

\(第二项=\int_{\Omega} \varepsilon(v): C: \varepsilon(u) d x = \int (BV)^T(DBU) dx = V^T(\int B^TDB dx) U = V^TKU\)

由于要求对任意v均成立，一定是\(M\ddot{U}+KU=0\)，这就是基于KM离散化的变形结果。在离散形式下可以做之前说的分时间步解线性方程组的模拟。

特征函数解法

为了避免前面说的每一步求矩阵逆，考虑把原PDE分解成多个ODE的“线性组合”的思路，每一步可以直接得到增量，不需要解方程。这个方法被称为Modal Analysis / Eigenfunction Expansion，这些ODE就是特征函数（模态），线性系数称为统本征频率。

具体来说，对于连续形式（强形式），找到一组函数满足：

1. \(u(x,t)\) = \(\Sigma \phi_i(x)q_i(t)\) 带回后满足原PDE
2. 这组函数满足正交性：\(\int \rho \phi_n\phi_m dx =0\)

第二个条件需要满足的原因是，这样带回PDE后，等式成立要求每一项对应相等，可以解耦成n个ODE问题。

关于为什么能找到，我查到是“线性、无阻尼、均匀材料或可正交化的系统”就能够做模态分离（x和t分离），所以第一步能找到；不能做的情况，就是不能用这个方法，不讨论。在这个前提下，得到的特征函数一定是正交的。因此流程就是设\(u(x,t)\) = \(\Sigma \phi_i(x)q_i(t)\)带回原PDE，然后分成n个ODE，加上边界条件，求解。

为什么保证正交我没看懂（），离散的形式好理解一些。离散下，\(M\ddot{U}+KU=0\)，设\(U(t)=\Phi q(t)\)，\(\Phi\)是特征矩阵，由n个列向量\(\phi_i\)组成。带回原式子：

\[M\Phi \ddot{q(t)} + K\Phi q(t) = 0 \]

想要拆出 \(K\phi_i q_i(t) = M\phi_i \ddot{q_i(t)}\)，首先求\(\Lambda\)使得 \(K\Phi = M \Phi \Lambda\)，这里\(\Phi\)是特征矩阵，\(\Lambda\)是对角阵（特征值），当K和M均对称正定，一定能找到，我们也只讨论KM均对称正定的情况。

假设\(\Lambda\)和\(\Phi\)对K归一化，使得\(\Phi^T M \Phi =I\)，则\(\Phi^T K \Phi = \Lambda\)，对原式子左乘\(\Phi^T\):

\[\Phi^T M\Phi \ddot{q(t)} + \Phi^T K\Phi q(t) = 0 \]

\[\ddot{q(t)} + \Lambda q(t) = 0 \]

由于\(\Lambda\)是对角阵，可以拆开成n个\(\ddot{q_i(t)} + \lambda_i q_i(t) = 0\)，是ODE问题。

神经特征函数

写不完了先去吃饭。