偏微分方程计算方法

2026-07-22

核心主题：有限差分法、基矢展开法（含 DVR）、B-样条、小波分析、有限元法、复杂边界条件处理

目录

1. 第一章：偏微分方程的有限差分法
2. 第二章：二维 Poisson 方程 — 五点差分与矩阵形式
3. 第三章：边界条件的三类处理
4. 第四章：傅里叶微分算符与 FFT
5. 第五章：基矢展开法求微分方程
6. 第六章：离散变量表示 DVR
7. 第七章：B-样条展开法
8. 第八章：小波分析与多小波展开
9. 第九章：有限元方法 FEM
10. 第十章：边界条件进阶 — 多极展开、周期性、Twisted
11. 第十一章：Jacobi 松弛迭代法
12. 总结与对比

第一章：偏微分方程的有限差分法

💡 物理直觉

偏微分方程（PDE）在物理中无处不在——Poisson 方程描述静电势、引力势；Schrödinger 方程描述量子系统的波函数；扩散方程描述热传导和物质输运。它们共同的特点是多变量 (x, y, z) 通过函数耦合在一起。

说人话：解一维问题就像走钢丝，解二维三维问题就像走钢丝网——你必须同时照顾好所有方向的"约束力"。

从一维到二维的推广

一维有限差分中，二阶导数离散化为：

$$ u''(x_i) \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{h^2} $$

在二维网格上，点 $(x_i, y_j)$ 的 Laplace 算子 $\nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$ 变成五点差分：

$$ \nabla^2 u\big|_{i,j} \approx \frac{u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4u_{i,j}}{h^2} $$

五点查分的具体推导：

二维 Laplace 算子是：

$$ \nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} $$ 在均匀网格上：

$$ x_i=x_0+ih,\qquad y_j=y_0+jh $$ 函数值记为：

$$ u_{i,j}=u(x_i,y_j) $$ 对 $x$ 方向二阶导：

$$ \frac{\partial^2 u}{\partial x^2}\bigg|_{i,j} \approx \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2} $$ 对 $y$ 方向二阶导：

$$ \frac{\partial^2 u}{\partial y^2}\bigg|_{i,j} \approx \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^2} $$ 两者相加：

$$ \nabla^2 u\big|_{i,j} \approx \frac{ u_{i+1,j} +u_{i-1,j} +u_{i,j+1} +u_{i,j-1} -4u_{i,j} }{h^2} $$ 这就是二维五点差分格式。

如果 x,y 方向上 网格间距不同

若：

$$ h_x\neq h_y $$ 则：

$$ \nabla^2 u\big|_{i,j} \approx \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h_x^2} + \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_y^2} $$ 整理：

$$ \nabla^2 u\big|_{i,j} \approx \frac{u_{i+1,j}+u_{i-1,j}}{h_x^2} + \frac{u_{i,j+1}+u_{i,j-1}}{h_y^2} - 2\left( \frac{1}{h_x^2} + \frac{1}{h_y^2} \right)u_{i,j} $$ 如果 $h_x=h_y=h$，就退化为五点公式。

📐 直观理解：一个点要"看"四个邻居

五点公式：

$$ \nabla^{2}u_{i,j} \approx \frac{u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j}}{h^{2}} $$ 可以写成：

$$ \nabla^{2}u_{i,j} \approx \frac{4}{h^2}\left(\frac{u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}}{4} -u_{i,j}\right) $$

也就是说 约等于 邻居的平均值 减去 自身 就是 $\nabla^{2}$

所以：

• 如果 $u_{i,j}$ 大于周围平均值，Laplace 值偏负；
• 如果 $u_{i,j}$ 小于周围平均值，Laplace 值偏正；
• 如果 $u_{i,j}$ 等于周围平均值，Laplace 值为零。

            u[i, j+1]
               ↑
u[i-1, j] ← u[i, j] → u[i+1, j]
               ↓
            u[i, j-1]

说人话：在格点上算二阶导数，就是拿上下左右四个邻居的值减去4倍自己，再除以格距平方。这本质上是在问："我和邻居们的平均值差了多少？"

以此 ， 可以解释 调和函数 可以导出的平均值性质：

$$ \nabla^{2}u = 0 \Longrightarrow u_{i,j} = \frac{u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}}{4} $$ 即，没有源的时候，一个点的值等于周围“邻居” 的平均值

同理 可以推广到 三维时候的 七点差分（就是看一共六个“邻居”）

第二章：二维 Poisson 方程 — 五点差分与矩阵形式

💡 物理直觉

泊松方程 $\nabla^2 u = -f$ 的物理意义：给定一个"源"分布 $f$（如电荷密度），求它产生的"势" $u$（如电势）。五点差分把这个连续问题转化为一个巨大的线性方程组。

$$ 4u_{i,j} -u_{i+1,j} -u_{i-1,j} -u_{i,j+1} -u_{i,j-1} = h^2 f_{i,j} $$

重要提醒：两种矩阵约定

这里有一个非常容易混乱的地方。

约定 A：矩阵里面保留 $1/h^2$

写成：

$$ A=\frac{1}{h^2}A_0 $$ 其中 $A_0$ 的主对角是 $4$，邻居位置是 $-1$。

此时右端应该是：

$$ b=f+\frac{\text{边界贡献}}{h^2} $$

约定 B：整体乘以 $h^2$

写成： $$ A_0\mathbf u=\mathbf b $$ 此时右端是： $$ b=h^2f+\text{边界贡献} $$ 你给的伪代码：

Python

b[k] = f[i][j] * h**2

b[k] += boundary_value

对应的是第二种约定，也就是矩阵中不带 $1/h^2$。

如果矩阵写成： $$ A=\frac1{h^2}A_0 $$ 那么 RHS 就不该写 $f h^2$，而应该写 $f$。

具体数值例子：3×3 内部格点

设求解区域内部有 3×3 = 9 个未知格点，边界值已知（Dirichlet 条件）。

步骤 1：把 2D 格点"拉直"成一维向量

设未知函数在格点 $(i,j)$ 处的值为 $u_{i,j}$。按行（或列）排列：

$$ \mathbf{u} = [u_{1,1}, u_{1,2}, u_{1,3}, u_{2,1}, u_{2,2}, u_{2,3}, u_{3,1}, u_{3,2}, u_{3,3}]^T $$

这就是投影转为一维的核心思想——多维问题的本质就是把高维指标压扁成一个长向量。

步骤 2：构建微分算符矩阵 A

$A = \frac{1}{h^2} A_{0}$

对内部 3×3 网格，离散 Laplace 算子对应的矩阵为 9×9 分块三对角形式。取 $h=1$ 简化说明：

$$ A = \frac{1}{h^2} \begin{bmatrix} 4 & -1 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ -1 & 4 & -1 & 0 & -1 & 0 & 0 & 0 & 0 \\ 0 & -1 & 4 & 0 & 0 & -1 & 0 & 0 & 0 \\ -1 & 0 & 0 & 4 & -1 & 0 & -1 & 0 & 0 \\ 0 & -1 & 0 & -1 & 4 & -1 & 0 & -1 & 0 \\ 0 & 0 & -1 & 0 & -1 & 4 & 0 & 0 & -1 \\ 0 & 0 & 0 & -1 & 0 & 0 & 4 & -1 & 0 \\ 0 & 0 & 0 & 0 & -1 & 0 & -1 & 4 & -1 \\ 0 & 0 & 0 & 0 & 0 & -1 & 0 & -1 & 4 \end{bmatrix} $$

说人话：矩阵每行最多 5 个非零元——主对角线是 4（自己），周围的 -1 对应上下左右四个邻居。矩阵是稀疏的、对称的、正定的。这种结构在工程中叫"块三对角"（block tridiagonal）。

步骤 3：边界条件进入右端向量

以左下角附近的内部点 $u_{1,1}$ 为例。

离散方程：

$$ 4u_{1,1} -u_{2,1} -u_{0,1} -u_{1,2} -u_{1,0} = h^2f_{1,1} $$ 其中：

• $u_{2,1}$、$u_{1,2}$ 是内部未知量；
• $u_{0,1}$、$u_{1,0}$ 是边界已知量。

把边界已知量移到右端：

$$ 4u_{1,1} -u_{2,1} -u_{1,2} = h^2f_{1,1} + u_{0,1} + u_{1,0} $$ 所以边界值是以正号加到 RHS 上的。

一般规律：

$$ b_{i,j} = h^2f_{i,j} + \sum_{\text{邻居为边界}} u_{\text{boundary}} $$ 这是在矩阵不带 $1/h^2$ 的约定下。

# 伪代码：构建右端向量并加入边界贡献
for i in range(1, Nx-1):
    for j in range(1, Ny-1):
        k = idx(i, j)  # 2D→1D 映射
        b[k] = f[i][j] * h**2  # 源项
        # 边界贡献
        if i == 1:     b[k] += u[0][j]       # 下边界
        if i == Nx-2:  b[k] += u[Nx-1][j]    # 上边界
        if j == 1:     b[k] += u[i][0]       # 左边界
        if j == Ny-2:  b[k] += u[i][Ny-1]    # 右边界

步骤 4：求解线性方程组 $$ A_{0}\mathbf{u} = \mathbf{b} $$

用 SciPy 稀疏求解器一行搞定：

from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve

# A 为稀疏矩阵，b 为右端向量
u_1d = spsolve(A, b)
u_2d = u_1d.reshape(Nx, Ny)  # 还原为二维

用Kronecker积写二维的Laplace矩阵！

二维矩阵其实可以用一维矩阵拼出来。

一维二阶差分矩阵： $$ T_N= \begin{bmatrix} 2 & -1 & 0 & \cdots & 0\\ -1 & 2 & -1 & \cdots & 0\\ 0 & -1 & 2 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & -1\\ 0 & 0 & 0 & -1 & 2 \end{bmatrix} $$ 对于二维 Poisson： $$ -\nabla^2 = -\partial_{xx} -\partial_{yy} $$ 离散矩阵可以写成： $$ A= \frac1{h^2} \left( I_y\otimes T_x + T_y\otimes I_x \right) $$ 如果 $N_x=N_y=3$，就会得到你写的 $9\times 9$ 块三对角矩阵。

这说明二维 Laplace 算子本质上是两个方向一维 Laplace 算子的张量和。

矩阵性质

对于 Dirichlet 边界条件，二维 Poisson 矩阵满足：

稀疏：每行最多 5 个非零元；

对称：

$$ A^T=A $$ 正定：

$$ \mathbf u^T A\mathbf u>0 $$ 只要 $\mathbf u\neq 0$。

正定性来自物理能量：

$$ \int |\nabla u|^2\,d\Omega>0 $$ 离散后对应：

$$ \mathbf u^TA\mathbf u \sim \sum_{\langle i,j\rangle}(u_i-u_j)^2 $$ 所以这个矩阵适合用共轭梯度法 CG

📊 矩阵规模与计算代价

⚠️ 关键认识：矩阵维数 = (格点数)²，二维问题的计算量随格点数增长极快。这就是为什么三维问题经常让人绝望——矩阵维度 = (格点数)³ 的平方！

💡 算法思路：高阶差分与精度权衡

可以采用第三章学过的拉格朗日插值法构造更高阶的差分公式（如5点、7点格式）。但在边界处只能退化为低阶格式，否则会用到界外的不存在的格点。

核心矛盾：格点间距越小精度越高，但二维三维很难取很小的间距（存储和计算都爆炸）。因此需要考虑：

• ✨ 对称性降维：利用宇称对称性可以减少一半计算空间
• ✨ 好量子数：量子力学中利用好量子数进行变量分离（如角动量投影）
• ⚠️ 一旦对称性破缺，计算量增加几个量级

第三章：边界条件的三类处理

PDE 的定解问题 = 控制方程 + 边界条件。

边界条件就是"物理系统在边界上被怎么约束"。不给定边界条件，解就不唯一——就像做菜不告诉我锅的温度，我咋知道该用多大火力？

否则解通常不唯一

三类边界条件对比

Neumann 边界条件的数值处理

Neumann 条件不像 Dirichlet 那样可以直接"削掉"边界格点——因为边界值本身是未知的。要引入额外的差分方程。

例如引入虚格点（ghost point）：在边界外对称放一个辅助点，利用中心差分 $\frac{u_{i+1} - u_{i-1}}{2h} = g$ 来消去虚点。

一维 Neumann 的 虚格点 推导

考虑一维 Poisson：

$$ -u''(x)=f(x),\qquad x\in[0,L] $$ 左边界给 Neumann：

$$ u'(0)=g $$ 在 $x=0$ 处引入虚点：

$$ x_{-1}=-h $$ 中心差分近似一阶导数：

$$ u'(0) \approx \frac{u_1-u_{-1}}{2h} $$ 令它等于 $g$：

$$ \frac{u_1-u_{-1}}{2h}=g $$ 解出虚点：

$$ u_{-1}=u_1-2hg $$ 现在在边界点 $x_0$ 上离散 PDE：

$$ -u''(0)=f_0 $$ 二阶导：

$$ u''(0) \approx \frac{u_1-2u_0+u_{-1}}{h^2} $$ 所以：

$$ -\frac{u_1-2u_0+u_{-1}}{h^2}=f_0 $$ 代入：

$$ u_{-1}=u_1-2hg $$ 得到：

$$ -\frac{u_1-2u_0+u_1-2hg}{h^2}=f_0 $$ 整理：

$$ -\frac{2u_1-2u_0-2hg}{h^2}=f_0 $$ 因此：

$$ \frac{2u_0-2u_1}{h^2} = f_0-\frac{2g}{h} $$ 这就是 Neumann 边界点对应的离散方程。

如果 $g=0$，即对称边界：

$$ u'(0)=0 $$ 则：

$$ u_{-1}=u_1 $$ 这意味着边界外的虚点是边界内点的镜像。

二维 Neumann 边界

假设左边界 $x=0$ 给：

$$ \frac{\partial u}{\partial x}(0,y_j)=g_j $$ 引入虚点 $u_{-1,j}$，有：

$$ \frac{u_{1,j}-u_{-1,j}}{2h}=g_j $$ 所以：

$$ u_{-1,j}=u_{1,j}-2hg_j $$ 如果边界点本身也是未知量 $u_{0,j}$，则在边界处的 Laplace 方程：

$$ -\nabla^2u=f $$ 写成：

$$ -\frac{ u_{1,j}+u_{-1,j}+u_{0,j+1}+u_{0,j-1}-4u_{0,j} }{h^2} = f_{0,j} $$ 代入虚点：

$$ u_{-1,j}=u_{1,j}-2hg_j $$ 得到：

$$ -\frac{ 2u_{1,j} -2hg_j +u_{0,j+1} +u_{0,j-1} -4u_{0,j} }{h^2} = f_{0,j} $$ 这就把虚点消掉了。

Robin 边界条件

Robin 条件：

$$ \alpha u+\beta\frac{\partial u}{\partial n}=g $$ 例如一维左边界：

$$ \alpha u_0+\beta u'(0)=g $$ 用中心差分：

$$ u'(0)\approx \frac{u_1-u_{-1}}{2h} $$ 代入：

$$ \alpha u_0 + \beta\frac{u_1-u_{-1}}{2h} = g $$ 解出虚点：

$$ u_{-1} = u_1-\frac{2h}{\beta}(g-\alpha u_0) $$ 即：

$$ u_{-1} = u_1-\frac{2hg}{\beta} + \frac{2h\alpha}{\beta}u_0 $$ 然后再代回 PDE 的差分方程。

纯Neumann问题的特殊性

如果整个边界都是 Neumann 条件，例如：

$$ -\nabla^2u=f \\ \space \\ \frac{\partial u}{\partial n} = g , on \space \partial \Omega $$ 则解只确定到一个常数。

因为如果 $u$ 是解，那么：

$$ u+C $$ 也是解。

所以矩阵会有零特征值，不是严格正定。

还需要兼容条件。

由散度定理：

$$ \int_\Omega -\nabla^2u\,d\Omega = -\int_{\partial\Omega} \frac{\partial u}{\partial n} \,dS $$ 因为：

$$ -\nabla^2u=f $$ 所以：

$$ \int_\Omega f\,d\Omega = -\int_{\partial\Omega}g\,dS $$ 如果这个条件不满足，纯 Neumann 问题无解。

如果满足，解存在但差一个常数，需要额外固定： $$ \int_\Omega u\,d\Omega=0 $$ 或指定某一点： $$ u(x_0)=0 $$

说人话：Dirichlet 就是"墙的温度给我定死了"，Neumann 就是"墙的热流给我定死了"。前者简单（少解几个方程），后者麻烦（多解几个方程）。

第四章：傅里叶微分算符与 FFT

💡 物理直觉

有限差分 就是 在实空间中用邻居差值近似导数

Fourier方法 是在 频率/动量空间 利用：

$$ \frac{d}{dx} e^{ikx} = ik e^{ikx} $$ 在格点空间做微分的"近邻相减"，在动量空间就是简单的乘法！

这是所有谱方法的精髓 （微分变成乘法）

离散 Fourier变换 DFT

在周期区间：

$$ x_j=jh,\qquad j=0,1,\dots,N-1 $$ 离散 Fourier 展开：

$$ u_j = \sum_{m=0}^{N-1} \tilde u_m e^{ik_m x_j} $$ 其中波数：

$$ k_m=\frac{2\pi}{L}m $$ 但是 DFT 的频率顺序通常是：

$$ m=0,1,\dots,\frac N2,-\frac N2+1,\dots,-1 $$ 在 NumPy 中：

Python
k = 2*np.pi*np.fft.fftfreq(N, d=h)

会自动给出正确频率。

在动量空间，一阶微分就是乘 $ik$，二阶微分就是乘 $-k^2$：

$$ \partial_x \leftrightarrow ik, \qquad \partial_x^2 \leftrightarrow -k^2 $$

再逆变换回来，就得到了格点上的微分近似。

Fourier 微分矩阵

令 DFT 矩阵为 $F$，则：

$$ \tilde{\mathbf u}=F\mathbf u $$ 反变换：

$$ \mathbf u=F^{-1}\tilde{\mathbf u} $$ 一阶导数在 Fourier 空间是乘 $ik$，所以：

$$ \widetilde{\mathbf u'} = \operatorname{diag}(ik_m)\tilde{\mathbf u} $$ 回到实空间：

$$ \mathbf u' = F^{-1}\operatorname{diag}(ik_m)F\mathbf u $$ 因此一阶微分矩阵是：

$$ D^{(1)} = F^{-1}\operatorname{diag}(ik_m)F $$ 二阶微分矩阵：

$$ D^{(2)} = F^{-1}\operatorname{diag}(-k_m^2)F $$ 因为：

$$ (ik_m)^2=-k_m^2 $$ 所以：

$$ D^{(2)} = D^{(1)}D^{(1)} $$ 注意：这个关系在 Fourier 谱方法里严格成立。

但普通有限差分中，中心差分的一阶矩阵平方，不一定等于常用的二阶差分矩阵。

一阶和二阶微分算符矩阵的关系

Fourier 方法的优雅之处：一阶矩阵 $D^{(1)}$ 和二阶矩阵 $D^{(2)}$ 是完全一致的——$D^{(2)} = D^{(1)}D^{(1)}$，即微分和积分互逆。

FFT 的计算效率

上一轮我们推导的傅里叶微分矩阵 $D=F^{-1}diag(ik)F$ ，在数学上完美无瑕，但如果直接按这个公式去编程，你的电脑会卡死。这段文字就是在告诉你：FFT 是专门来拯救这个困局的“救世主”

FFT发现了 一个极其简单但有效的规律 ： ==偶数项和奇数项可以分开算==

• 一个 NN 点的大变换，可以拆成 两个 N/2点的小变换（一个只负责偶数位置，一个只负责奇数位置）。
• 而这两个 N/2N/2 点的小变换，又可以继续拆成 四个 N/4 点的更小变换……
• 直到拆成单个点（单个点的傅里叶变换就是它自己）。

这就是“分治（Divide and Conquer）”。 通过这种层层拆分，原本需要把每个格点和每个频率都配对的 N2N2 次运算，变成了只需在每一层（共 log⁡NlogN 层）把结果合并起来（合并时只涉及少量运算）的 Nlog⁡NNlogN 次运算。

说人话：FFT 把"每个格点和每个动量分量两两配对"的暴力算法，变成了精巧的分治算法。原来 N² 的事现在 N log N 就办完了——当 N 很大时，这个差距是天壤之别。

理论上的 Fourier微分矩阵 虽然是稠密的（计算量巨大）但实际我们从来不显式地去构建这个矩阵。我们直接用FFT代替矩阵乘法，把计算复杂度从“平方级”降为“线性对数级”，从而使得谱方法能够处理 百万级网格的工程问题

🔬 二维 Fourier 微分

二维周期函数：

$$ u(x,y) = \sum_{k_x,k_y} \tilde u_{k_x,k_y} e^{i(k_xx+k_yy)} $$ 则：

$$ \nabla^2u = \frac{\partial^2u}{\partial x^2} + \frac{\partial^2u}{\partial y^2} $$ 在 Fourier 空间变成：

$$ \widehat{\nabla^2u} = -(k_x^2+k_y^2)\tilde u $$ 所以算法是：

1. 对 $u$ 做二维 FFT；
2. 每个 Fourier 分量乘上 $-(k_x^2+k_y^2)$；
3. 做逆 FFT。

二维情况使用张量积形式：依次对每行、每列做一维变换。

import numpy as np

def fourier_laplacian_2d(u, hx, hy=None):
    """
    用 FFT 计算二维 Laplace 算子。
    要求周期边界条件。
    """
    if hy is None:
        hy = hx

    Nx, Ny = u.shape

    kx = 2*np.pi*np.fft.fftfreq(Nx, d=hx)
    ky = 2*np.pi*np.fft.fftfreq(Ny, d=hy)

    KX, KY = np.meshgrid(kx, ky, indexing='ij')

    multiplier = -(KX**2 + KY**2)

    u_hat = np.fft.fft2(u)
    lap_hat = multiplier * u_hat

    return np.fft.ifft2(lap_hat).real

🎯 FFT 的应用拓展

Fourier方法 默认周期边界条件！！！

也就是说：

$$ u(0)=u(L) $$ 如果函数不周期，强行用 Fourier 会导致边界跳跃，从而产生 Gibbs 现象。

对于非周期边界，可以用：

• sine transform：适合 Dirichlet；
• cosine transform：适合 Neumann；
• Chebyshev 谱方法：适合有限区间非周期问题。

具体应用：

• 图像处理：二维 FFT 是 JPEG 压缩、卷积神经网络的基础
• 信号处理：去噪（低通滤波）、边缘检测（高通滤波）
• 零填充（zero-padding）：平滑处理
• 共振态展宽：$A(t) = A e^{i\omega t} e^{-\Gamma t/2}$ → 傅里叶变换得到 Lorentz 线形，展宽 $\Gamma$ 对应寿命 $\tau = 1/\Gamma$

用 FFT 解 Poisson 方程

考虑周期边界下：

$$ -\nabla^2u=f $$ Fourier 空间：

$$ -\widehat{\nabla^2u} = (k_x^2+k_y^2)\hat u = \hat f $$ 所以：

$$ \hat u = \frac{\hat f}{k_x^2+k_y^2} $$ 但当：

$$ k_x=k_y=0 $$ 分母为零。

这是因为周期 Poisson 方程的解只确定到一个常数。

必须要求：

$$ \int f\,d\Omega=0 $$ 也就是：

$$ \hat f_{0,0}=0 $$ 然后可以人为设定：

$$ \hat u_{0,0}=0 $$ 即选择零平均势。

FFT 与共振展宽

如果一个信号：

$$ A(t)=A_0e^{-i\omega_0t}e^{-\Gamma t/2} $$ 它的 Fourier 变换近似为：

$$ \tilde A(\omega) = \int_0^\infty A_0e^{-i\omega_0t}e^{-\Gamma t/2}e^{i\omega t}\,dt $$ 整理指数：

$$ \tilde A(\omega) = A_0 \int_0^\infty e^{-\left(\frac{\Gamma}{2}+i(\omega_0-\omega)\right)t} \,dt $$ 积分：

$$ \tilde A(\omega) = \frac{A_0} { \frac{\Gamma}{2}+i(\omega_0-\omega) } $$ 谱强度：

$$ |\tilde A(\omega)|^2 = \frac{|A_0|^2} { (\omega-\omega_0)^2+\left(\frac{\Gamma}{2}\right)^2 } $$ 这就是 Lorentz 线形。

$$ \Gamma $$ 越大，谱线越宽，寿命越短：

$$ \tau\sim \frac1{\Gamma} $$ 若带 $\hbar$，则：

$$ \tau\sim \frac{\hbar}{\Gamma} $$

第五章：基矢展开法求微分方程

💡 物理直觉

量子力学课上学过：波函数可以在任意一组完备基下展开。计算物理中我们就用这个思路——不是直接在格点上离散，而是把解写成一组合适基函数的线性组合，然后把 PDE 转化为矩阵本征值问题。

基本公式推导

将波函数展开为基矢 $\phi_j$ 的线性组合：

$$ \psi_i = \sum_j c_{ij} \phi_j $$

代入本征方程 $H\psi = E\psi$：

$$ \sum_m H \phi_m \, c_{im} = E \sum_j c_{ij} \phi_j $$

左乘 $\phi_k^*$ 并积分（投影）：

$$ \sum_m \langle \phi_k | H | \phi_m \rangle \, c_{im} = E \, c_{ik} $$

写成矩阵形式：

$$ \begin{bmatrix} H_{11} & H_{12} & \cdots & H_{1N} \\ H_{21} & H_{22} & \cdots & H_{2N} \\ \vdots & \vdots & \ddots & \vdots \\ H_{N1} & H_{N2} & \cdots & H_{NN} \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_N \end{bmatrix} = E \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_N \end{bmatrix} $$

说人话：本来是一个微分方程，经过"选定基函数→展开→投影→积分"这一套，变成了一个矩阵对角化问题。这就是Ritz 变分法的数值实现。

对于基态能量，有限基截断给出的能量通常是上界：

$$ E_0^{(N)}\ge E_0 $$ 随着基函数数目增大：

$$ E_0^{(N)}\to E_0 $$

实际应用：三维谐振子基矢

3D 笛卡尔谐振子基

三维笛卡尔谐振子基：

$$ \phi_{n_xn_yn_z}(x,y,z) = \phi_{n_x}(x) \phi_{n_y}(y) \phi_{n_z}(z) $$ 其中一维谐振子基函数：

$$ \phi_n(x) = N_n H_n(\xi)e^{-\xi^2/2} $$ $H_n$ 是 Hermite 多项式，$b$ 是振子长度参数。

总壳层量子数：

$$ N=n_x+n_y+n_z $$ 如果截断到：

$$ N\le N_{\max} $$ 则基函数数目是：

$$ \sum_{N=0}^{N_{\max}} \frac{(N+1)(N+2)}{2} = \frac{(N_{\max}+1)(N_{\max}+2)(N_{\max}+3)}{6} $$ 这个数目随着 $N_{\max}$ 迅速增长。

轴对称体系（柱坐标）

对于轴对称系统，通常用柱坐标：

$$ (\rho,\varphi,z) $$ 基函数可以分离为：

$$ \Phi(\rho,\varphi,z,\sigma) = R_{n_r}^{|m|}(\rho) e^{im\varphi} Z_{n_z}(z) \chi_\sigma $$ 其中：

• $Z_{n_z}(z)$：Hermite 基；
• $R_{n_r}^{|m|}(\rho)$：associated Laguerre 基；
• $e^{im\varphi}$：角向部分；
• $\chi_\sigma$：自旋部分。

常用截断：

$$ N=2n_r+n_z+|m| $$ 保留：

$$ N\le N_{\max} $$

好量子数为什么能降维？

如果 Hamiltonian 和某个算符对易： $$ [H,\hat Q]=0 $$ 那么 $\hat Q$ 的本征值就是好量子数。

例如轴对称情况下，角动量 $z$ 分量 $m$ 可能守恒。

于是不同 $m$ 的态不会耦合。

Hamiltonian 矩阵可以分块： $$ H= \begin{bmatrix} H_{m=0} & 0 & 0\\ 0 & H_{m=1} & 0\\ 0 & 0 & H_{m=2} \end{bmatrix} $$ 原来需要对角化一个巨大矩阵，现在可以分别对角化很多小矩阵。

这就是利用对称性降维的核心。

一旦对称性破缺，原本的块之间开始耦合，矩阵规模会暴涨。

💡 关键认识

• 基函数自带边界条件：不像有限差分需要显式处理边界，基函的选择就已隐含了边界行为（如谐振子基在无穷远处指数衰减）
• 微分可解析计算：谐振子基的微分有解析表达式，不需要数值近似
• 截断误差：基函数个数 N 有限，所以是近似方法。N 越大越精确
• 适用场景：形状规则的束缚态问题效果极好；复杂形状或连续谱则不够精确

流程对比

都是把无限维问题 最后 变成 有限维问题

有限差分法：  格点离散 → 差分近似微分 → 矩阵方程求解 （在格点上表示函数）
基矢展开法：  选基函数 → 计算矩阵元 → 矩阵对角化 （用系数表示函数）

第六章：离散变量表示 DVR

💡 物理直觉

DVR（Discrete Variable Representation）是一种同时兼具基空间和格点空间优点的表示方法。它让你"吃着碗里的看着锅里的"——基函数保证了完备性，而格点表示让矩阵元计算极其简单。

==核心：动能像基空间一样可以精确/解析处理，势能像格点空间一样近似对角。==

DVR 的核心思想

给定一组正交基函数 $\phi_n(x)$ 和高斯积分点（gauge points）$\{x_\alpha\}$：

1. 用这组基构造投影算符 P
2. 在格点 $\{x_\alpha\}$ 上定义 DVR 基函数：$f_\alpha(x)$，满足 $f_\alpha(x_\beta) = \delta_{\alpha\beta}$（在自家格点为 1，别家格点为 0）
3. 关键性质：

• 动能矩阵元：有解析表达式，计算简单
• 势能矩阵元：对角！$V_{\alpha\beta} \approx V(x_\alpha)\delta_{\alpha\beta}$

例子

均匀格点上常用 sinc-DVR。

假设区间长度为 $L$，格点数为 $N$，格距：

$$ \Delta x=\frac{L}{N} $$ 对于无限 sinc-DVR，常见二阶导数矩阵是：

$$ \left(\frac{d^2}{dx^2}\right)_{ij} = \begin{cases} -\dfrac{\pi^2}{3\Delta x^2}, & i=j\\[8pt] -\dfrac{2(-1)^{i-j}}{\Delta x^2(i-j)^2}, & i\neq j \end{cases} $$ 于是动能：

$$ T=-\frac{\hbar^2}{2m}\frac{d^2}{dx^2} $$ 对应：

$$ T_{ij} = \frac{\hbar^2}{2m} \begin{cases} \dfrac{\pi^2}{3\Delta x^2}, & i=j\\[8pt] \dfrac{2(-1)^{i-j}}{\Delta x^2(i-j)^2}, & i\neq j \end{cases} $$ 周期 sinc-DVR 或 Fourier-DVR 的公式会略有不同，比如会出现：

$$ \sin^2\left(\frac{\pi(i-j)}{N}\right) $$ 这正是你代码里的结构。

需要注意：

不同文献里 $T$、$-d^2/dx^2$、$-\frac12d^2/dx^2$ 的定义可能不同，所以符号和系数要仔细核对。

# DVR 动能矩阵元的典型计算框架
# 以 sinc-DVR 为例（均匀格点）
def dvr_kinetic_1d(N, L):
    """一维 sinc-DVR 动能矩阵，区间 [-L/2, L/2]"""
    dx = L / N
    T = np.zeros((N, N))
    for i in range(N):
        for j in range(N):
            if i == j:
                T[i, j] = np.pi**2 / (6 * dx**2) * (1 + 2/N**2)
            else:
                d = i - j
                T[i, j] = (-1)**d / (dx**2 * np.sin(np.pi * d / N)**2)
    return T

从 FBR（有限基表示） 到 DVR

给定正交基函数：

$$ \phi_n(x),\qquad n=0,1,\dots,N-1 $$ 定义投影算符：

$$ P=\sum_{n=0}^{N-1}|\phi_n\rangle\langle \phi_n| $$ 这个投影算符的意思是： 我们只在这 $N$ 个基函数张成的有限维空间里面近似真实波函数。

选择高斯积分点和权重

选一组积分点：

$$ x_\alpha,\qquad \alpha=0,1,\dots,N-1 $$ 以及对应权重：

$$ w_\alpha $$ 满足高斯积分近似：

$$ \int f(x)\,dx \approx \sum_{\alpha=0}^{N-1} w_\alpha f(x_\alpha) $$ 你讲义里写的是 gauge points，更准确地说一般叫：

$$ \text{quadrature points} $$ 或者：

$$ \text{DVR grid points} $$ 如果这组点和权重选得好，那么对于基函数乘积：

$$ \int \phi_m^*(x)\phi_n(x)\,dx $$ 可以用求和精确或高精度表示：

$$ \int \phi_m^*(x)\phi_n(x)\,dx \approx \sum_{\alpha=0}^{N-1} w_\alpha \phi_m^*(x_\alpha)\phi_n(x_\alpha) = \delta_{mn} $$

DVR 变换矩阵

定义矩阵：

$$ U_{n\alpha}=\sqrt{w_\alpha}\phi_n(x_\alpha) $$ 这里：

• $n$：基函数指标；
• $\alpha$：格点指标。

因为：

$$ \sum_\alpha w_\alpha \phi_m^*(x_\alpha)\phi_n(x_\alpha) \approx \delta_{mn} $$ 所以：

$$ \sum_\alpha U_{m\alpha}^*U_{n\alpha} = \delta_{mn} $$ 也就是：

$$ UU^\dagger \approx I $$ 在理想情况下，$U$ 是酉矩阵或正交矩阵。

DVR 基函数的构造

DVR 基函数定义为：

$$ |\chi_\alpha\rangle = \sum_{n=0}^{N-1} U_{n\alpha}^* |\phi_n\rangle $$ 如果基函数都是实函数，可以省略复共轭：

$$ |\chi_\alpha\rangle = \sum_{n=0}^{N-1} U_{n\alpha} |\phi_n\rangle $$

写成坐标形式：

$$ \chi_\alpha(x) = \sum_{n=0}^{N-1} \sqrt{w_\alpha}\phi_n(x_\alpha)\phi_n(x) $$ 这个 $\chi_\alpha(x)$ 就是 DVR 基函数。

DVR 基函数为什么具有“插值性”？

我们希望它满足类似：

$$ f_\alpha(x_\beta)=\delta_{\alpha\beta} $$ 也就是：

• 在自己的格点上为 1；
• 在别人的格点上为 0。

来推导。

在格点 $x_\beta$ 处：

$$ \chi_\alpha(x_\beta) = \sum_n \sqrt{w_\alpha} \phi_n(x_\alpha) \phi_n(x_\beta) $$ 利用矩阵 $U$ 的正交性：

$$ \sum_n U_{n\alpha}U_{n\beta} = \delta_{\alpha\beta} $$ 代入：

$$ U_{n\alpha} = \sqrt{w_\alpha}\phi_n(x_\alpha) $$ 所以：

$$ \sum_n \sqrt{w_\alpha}\phi_n(x_\alpha) \sqrt{w_\beta}\phi_n(x_\beta) = \delta_{\alpha\beta} $$ 两边除以 $\sqrt{w_\beta}$：

$$ \sum_n \sqrt{w_\alpha}\phi_n(x_\alpha) \phi_n(x_\beta) = \frac{\delta_{\alpha\beta}}{\sqrt{w_\beta}} $$ 左边正是 $\chi_\alpha(x_\beta)$。

因此：

$$ \chi_\alpha(x_\beta) = \frac{\delta_{\alpha\beta}}{\sqrt{w_\beta}} $$ 如果重新定义归一化后的 DVR cardinal 函数：

$$ f_\alpha(x)=\sqrt{w_\alpha}\chi_\alpha(x) $$ 那么：

$$ f_\alpha(x_\beta) = \sqrt{w_\alpha} \frac{\delta_{\alpha\beta}}{\sqrt{w_\beta}} $$ 当 $\alpha=\beta$ 时：

$$ f_\alpha(x_\alpha)=1 $$ 当 $\alpha\neq\beta$ 时：

$$ f_\alpha(x_\beta)=0 $$ 所以：

$$ \boxed{ f_\alpha(x_\beta)=\delta_{\alpha\beta} } $$ 这就是 DVR 的核心插值性质。

势能矩阵为什么近似对角？

考虑势能矩阵：

$$ V_{\alpha\beta}^{\text{DVR}} = \langle \chi_\alpha|V|\chi_\beta\rangle $$ 写成积分：

$$ V_{\alpha\beta} = \int \chi_\alpha^*(x)V(x)\chi_\beta(x)\,dx $$ 用高斯积分近似：

$$ V_{\alpha\beta} \approx \sum_\gamma w_\gamma \chi_\alpha^*(x_\gamma) V(x_\gamma) \chi_\beta(x_\gamma) $$ 利用：

$$ \chi_\alpha(x_\gamma) = \frac{\delta_{\alpha\gamma}}{\sqrt{w_\gamma}} $$ 以及：

$$ \chi_\beta(x_\gamma) = \frac{\delta_{\beta\gamma}}{\sqrt{w_\gamma}} $$ 代入：

$$ V_{\alpha\beta} \approx \sum_\gamma w_\gamma \left( \frac{\delta_{\alpha\gamma}}{\sqrt{w_\gamma}} \right) V(x_\gamma) \left( \frac{\delta_{\beta\gamma}}{\sqrt{w_\gamma}} \right) $$ 整理：

$$ V_{\alpha\beta} \approx \sum_\gamma \delta_{\alpha\gamma}\delta_{\beta\gamma} V(x_\gamma) $$ 只有当：

$$ \alpha=\beta=\gamma $$ 时有贡献。

所以：

$$ \boxed{ V_{\alpha\beta}^{\text{DVR}} \approx V(x_\alpha)\delta_{\alpha\beta} } $$ 这就是 DVR 最重要的结果。

势能矩阵不需要做复杂积分，直接取格点上的势能值。

动能矩阵如何处理？

动能算符为：

$$ T=-\frac{\hbar^2}{2m}\frac{d^2}{dx^2} $$ 在 FBR 中，动能矩阵是：

$$ T_{mn}^{\text{FBR}} = \left\langle \phi_m \left| -\frac{\hbar^2}{2m}\frac{d^2}{dx^2} \right| \phi_n \right\rangle $$ 这通常可以解析算出，或者用高精度积分算出。

然后用变换矩阵 $U$ 变到 DVR：

$$ T^{\text{DVR}} = U^\dagger T^{\text{FBR}}U $$ 即：

$$ T_{\alpha\beta}^{\text{DVR}} = \sum_{m,n} U_{m\alpha}^* T_{mn}^{\text{FBR}} U_{n\beta} $$ 所以 Hamiltonian 矩阵变成：

$$ H_{\alpha\beta}^{\text{DVR}} = T_{\alpha\beta}^{\text{DVR}} + V(x_\alpha)\delta_{\alpha\beta} $$ 最终求解：

$$ \boxed{ H\mathbf \psi=E\mathbf \psi } $$

DVR 的优势

📚 代表性工作

A. Bulgac and M.M. Forbes, Phys. Rev. C 87, 051301(R) (2013) — DVR 方法在核物理中的应用典范。

第七章：B-样条展开法

💡 物理直觉

有限差分用均匀格点，但有些地方（如原子核附近）波函数剧烈变化，有些地方（远离核）变化平缓。为什么不能"哪里变化快就多放格点，哪里变化慢就少放格点"？B-样条就是用来干这个的——它支持非均匀格点。

B-样条的定义

M 阶 B-样条 $B_i^M(x)$ 是在节点序列 $\{x_i, x_{i+1}, \dots, x_{i+M}\}$ 上由递推关系构造的分段多项式。在计算物理中常用的性质：

• 紧支撑：只在 M 个区间上非零
• 构成完备基
• 微分有解析递推公式
• 节点可以任意分布（边界加密！）

$M$：B-样条的阶数；

多项式次数为：

$$ p=M-1 $$ 例如：

• $M=1$：分段常数；
• $M=2$：分段线性；
• $M=3$：分段二次；
• $M=4$：分段三次 cubic spline。

一阶 B-样条定义为：

$$ B_i^1(x) = \begin{cases} 1, & t_i\le x

$$ B_i^M(x) = \frac{x-t_i}{t_{i+M-1}-t_i}B_i^{M-1}(x) + \frac{t_{i+M}-x}{t_{i+M}-t_{i+1}}B_{i+1}^{M-1}(x) $$ 如果分母为零，则对应项定义为零。

由递推关系可以看出：

$$ B_i^M(x) $$ 只在区间：

$$ [t_i,t_{i+M}] $$ 上非零。

所以它是紧支撑的。

紧支撑意味着矩阵是稀疏/带状的。

如果两个 B-样条支撑区间不重叠，那么：

$$ \int B_i^M(x)B_j^M(x)\,dx=0 $$ 因此重叠矩阵和 Hamiltonian 矩阵都是带状矩阵。

用 B-样条求解一维 Schrödinger 方程

考虑一维定态 Schrödinger 方程：

$$ \left[ -\frac{\hbar^2}{2m}\frac{d^2}{dx^2} + V(x) \right]\psi(x) = E\psi(x) $$ 将波函数展开：

$$ \psi(x) = \sum_j c_j B_j^M(x) $$ 代入方程：

$$ \left[ -\frac{\hbar^2}{2m}\frac{d^2}{dx^2} + V(x) \right] \sum_j c_jB_j^M(x) = E\sum_j c_jB_j^M(x) $$ 利用线性性：

$$ \sum_j c_j \left[ -\frac{\hbar^2}{2m}\frac{d^2B_j^M}{dx^2} + V(x)B_j^M(x) \right] = E\sum_jc_jB_j^M(x) $$ 左乘测试函数 $B_i^M(x)$，并积分：

$$ \int B_i^M(x) \sum_j c_j \left[ -\frac{\hbar^2}{2m}\frac{d^2B_j^M}{dx^2} + V(x)B_j^M(x) \right]dx = E \int B_i^M(x) \sum_jc_jB_j^M(x)dx $$ 交换求和与积分：

$$ \sum_jc_j \int B_i^M \left[ -\frac{\hbar^2}{2m}\frac{d^2B_j^M}{dx^2} + V B_j^M \right]dx = E\sum_jc_j \int B_i^M B_j^M dx $$ 定义：

$$ H_{ij} = \int B_i^M \left[ -\frac{\hbar^2}{2m}\frac{d^2}{dx^2} + V \right] B_j^M dx $$ 以及重叠矩阵：

$$ S_{ij}= \int B_i^M(x)B_j^M(x)\,dx $$ 于是：

$$ \boxed{ H\mathbf c=ES\mathbf c } $$ 这是广义本征值问题。

动能矩阵的更稳定写法

直接算：

$$ \int B_i\left(-\frac{d^2B_j}{dx^2}\right)dx $$ 有时候不如分部积分稳定。

考虑动能项：

$$ T_{ij} = -\frac{\hbar^2}{2m} \int B_i(x)B_j''(x)\,dx $$ 分部积分：

$$ \int B_i B_j''dx = \left[B_iB_j'\right]_{\text{boundary}} - \int B_i'B_j'dx $$ 如果 Dirichlet 边界下基函数或波函数边界为零，边界项消失：

$$ \left[B_iB_j'\right]_{\text{boundary}}=0 $$ 所以：

$$ T_{ij} = \frac{\hbar^2}{2m} \int B_i'(x)B_j'(x)\,dx $$ 因此：

$$ \boxed{ H_{ij} = \frac{\hbar^2}{2m}\int B_i'B_j'dx + \int B_i V B_j dx } $$ 这是实际计算中常用的弱形式。

为什么出现重叠矩阵 $S$？

普通正交基满足：

$$ \int \phi_i\phi_jdx=\delta_{ij} $$ 所以本征方程是：

$$ H\mathbf c=E\mathbf c $$ 但 B-样条通常不是正交基：

$$ \int B_iB_jdx\neq \delta_{ij} $$ 所以右边出现：

$$ S\mathbf c $$ 即：

$$ H\mathbf c=ES\mathbf c $$

不建议显式求 $S^{-1}$

讲义里写：

$$ (S^{-1}H)\mathbf c=E\mathbf c $$ 数学上可以这样写，但数值上不推荐显式求逆。

因为：

1. $S^{-1}H$ 可能破坏对称性；
2. 显式求逆误差大；
3. 矩阵可能变稠密。

更好的做法是直接解广义本征值问题：

Python

from scipy.linalg import eigh

eigvals, eigvecs = eigh(H, S)

如果一定要转成标准形式，可以做 Cholesky 分解：

$$ S=LL^T $$ 令：

$$ \mathbf c=L^{-T}\mathbf y $$ 代入：

$$ H L^{-T}\mathbf y = E S L^{-T}\mathbf y $$ 因为：

$$ S=L L^T $$ 所以：

$$ S L^{-T}=L $$ 左乘 $L^{-1}$：

$$ L^{-1}H L^{-T}\mathbf y=E\mathbf y $$ 得到标准对称本征值问题：

$$ \boxed{ \tilde H\mathbf y=E\mathbf y } $$ 其中：

$$ \tilde H=L^{-1}HL^{-T} $$

转换为标准本征值问题需要求 $B^{-1}$（最直观）

$$ (B^{-1} H) \mathbf{c} = E \mathbf{c} $$

说人话：B-样条和基矢展开法类似，但基函数不是解析多项式，而是分段拼接的函数段。矩阵方程中多了一个重叠矩阵 B，需要求逆转换。

B-样条节点怎么选？

常见节点分布：

1. 均匀节点

$$ t_i=t_0+ih $$ 适合变化均匀的问题。

2. 指数加密节点

原点附近密，远处稀：

$$ t_i=t_{\min} + \frac{e^{ai}-1}{e^{aN}-1} (t_{\max}-t_{\min}) $$ 适合原子、核、束缚态径向问题。

3. 分段节点

某些区域加密，其他区域粗化。

📊 精度与阶数的关系

随着 B-样条阶数 M 增大，精度单调改善。但需要注意 Nyquist 准则：格点间距不能太大，否则无法分辨波函数的振荡细节。

# 伪代码：B-样条展开法求解一维 Schrödinger 方程
# 1. 构造节点序列（边界加密）
knots = make_knots(x_min, x_max, n_basis, order=M)
# 2. 计算重叠矩阵 B 和 Hamiltonian 矩阵 H
B = overlap_matrix(knots, M)
H = hamiltonian_matrix(knots, M, potential_func)
# 3. 广义本征值问题 → 标准本征值问题
from scipy.linalg import eigh
B_inv = np.linalg.inv(B)
eigvals, eigvecs = eigh(B_inv @ H)

🔬 应用实例：形变体系的能级分布

在核物理中，B-样条被广泛用于求解轴对称形变势阱中的单粒子能级。利用对称性（角动量投影）可以极大降低矩阵维度。思考：对称性破缺后好量子数如何变化？

第八章：小波分析与多小波展开

💡 物理直觉

小波分析是一种多分辨率表示方法——对于波函数中振荡剧烈的区域，自动用精细的基来描述；对于平缓的区域，用粗糙的基来描述。这就像地图：城市用1:10000，沙漠用1:100000。

核心

Fourier 基函数是全局的：

$$ e^{ikx} $$ 它知道频率，但不知道这个频率出现在空间哪里。

小波基函数是局域的：

$$ \psi_{jk}(x) = 2^{j/2}\psi(2^jx-k) $$ 它同时有：

• 空间位置；
• 分辨率尺度。

其中：

• $j$：尺度，也叫 level；
• $k$：位置；
• $2^j$：缩放因子。
• $2^{j/2}$ 归一化因子

小波的两种表示方式

小波基具有两个核心操作：

• 平移 (Translation)：移动基函数的中心位置
• 伸缩 (Dilation)：缩放基函数的宽度

基函数族：$\psi_{jk}(x) = 2^{j/2} \psi(2^j x - k)$

其中 j 控制分辨率级别（level），k 控制中心位置。

多小波展开的核心思想

在 n-level 空间上展开波函数：

$$ f(x) = \sum_i s_i^n \phi_i^n(x) + \sum_{l=0}^{n-1} \sum_i d_i^l \psi_i^l(x) $$

其中：

• $s_i^n$ 是 n-level 的平滑系数（大尺度信息）
• $d_i^l$ 是各 level 的细节系数（逐级精细信息）
• 系数较小的 $d_i^l$ 可以截断 → 非均匀格点自动加密

核心：大尺度轮廓 + 各层细节修正

多小波是什么？

普通小波每个格子上通常一个 scaling function。

多小波 multiwavelet 在每个区间上使用多个基函数，比如 Legendre 多项式：

$$ \phi_0(x),\phi_1(x),\dots,\phi_{p-1}(x) $$ 这样每个单元内部不是只用一个常数，而是用高阶多项式。

优点：

• 高阶精度；
• 局部支撑；
• 自适应；
• 正交性好；
• 适合并行。

系数如何计算？

如果基函数正交，则：

$$ s_i^0=\langle \phi_i^0|f\rangle = \int \phi_i^0(x)f(x)\,dx $$ 如果某些细节系数很小：

$$ |d_i^l|<\epsilon $$ 就可以丢掉。

这就是自适应压缩。

🚀 算法创新

如果某个区域函数平滑，那么高层小波系数很小。

因为高层小波描述局部细节，而平滑区域没有多少局部细节。

所以： $$ d_i^l\approx 0 $$ 这些系数可以删掉。

如果某个区域函数剧烈变化，高层小波系数大，必须保留。

因此：

小波展开自动在变化快的地方保留更多自由度，在变化慢的地方保留更少自由度。

小波方法的核心创新：将微分方程转化为积分方程，迭代求解。这是代表性工作 MADNESS —— green方法啥的（G.I. Fann 等，ORNL）的核心思想。

说人话：小波方法=图像压缩算法的数学推广。JPEG 压缩扔掉高频细节，小波方法也扔掉小的展开系数，结果就是"波函数振荡的地方格点自动加密，平缓的地方格点自动稀疏"——无需人工干预！

流程类比

┌─────────────────────────────────────────────┐
│  输入: 微分方程 + 边界条件                     │
│    ↓                                        │
│  转换为积分方程                               │
│    ↓                                        │
│  多小波基展开 (n-level)                       │
│    ↓                                        │
│  截断小系数 → 自适应非均匀格点                  │
│    ↓                                        │
│  迭代求解 (如 GMRES)                          │
│    ↓                                        │
│  输出: 数值解 + 误差估计                       │
└─────────────────────────────────────────────┘

高性能实现

• TBB（Threading Building Blocks）：多线程并行
• Elemental：全局矩阵的分块分布式计算

第九章：有限元方法 FEM

💡 物理直觉

有限差分法在处理不规则形状边界时非常痛苦——网格线总不能完美贴合复杂边界。有限元法（FEM）的解决方案："把复杂蛋糕切成简单小块，在每一块上做简单近似，最后拼起来。"

有限差分法的想法是：

在规则网格上，用上下左右差分近似导数。

但 FEM 的想法完全不同：

不一定用规则网格，而是把区域切成很多小块。
每个小块里用简单函数近似 $u$。
最后把所有小块拼成一个整体。

FEM 不直接用差分近似二阶导数。

FEM 会说：

我不知道真实解 $u(x)$，但我可以用一堆简单函数拼出来。

令：

$$ u_h(x)=\sum_j U_j N_j(x) $$ 其中：

• $u_h(x)$：数值近似解；
• $U_j$：第 $j$ 个节点的未知值；
• $N_j(x)$：第 $j$ 个节点对应的基函数，也叫形函数 shape function。

它满足：

$$ N_1(x_1)=1 $$ 也就是说：

$N_1$ 在自己的节点为 1，在相邻节点为 0。

类似地：

$$ N_i(x_j)=\delta_{ij} $$ 这跟 DVR 里面的 cardinal property 有点像。

原方程是：

$$ -u''=f $$ 如果直接代入：

$$ u_h=\sum_j U_jN_j $$ 会出现二阶导数 $N_j''$。

但一维线性帽子函数是分段线性的，它的一阶导数是分段常数，二阶导数在节点处会有奇异性。

所以 FEM 不喜欢直接处理二阶导数。

于是它使用一个技巧：

把方程乘以测试函数，再积分，通过分部积分把二阶导数降成一阶导数。

这就叫 弱形式 weak form。

原方程：

$$ -u''(x)=f(x) $$ 取一个测试函数 $v(x)$，两边乘以 $v(x)$：

$$ -u''(x)v(x)=f(x)v(x) $$ 在整个区间积分：

$$ \int_0^L -u''(x)v(x)\,dx = \int_0^L f(x)v(x)\,dx $$ 所以：

$$ -\int_0^L u''v\,dx = -\left[u'v\right]_0^L + \int_0^L u'v'\,dx $$ 因此：

$$ \int_0^L u'v'\,dx - \left[u'v\right]_0^L = \int_0^L fv\,dx $$ 如果是 Dirichlet 边界：

$$ u(0)=0,\qquad u(L)=0 $$ 测试函数也取边界为零：

$$ v(0)=0,\qquad v(L)=0 $$ 于是边界项：

$$ \left[u'v\right]_0^L=0 $$ 所以得到：

$$ \boxed{ \int_0^L u'(x)v'(x)\,dx = \int_0^L f(x)v(x)\,dx } $$ 这就是弱形式。

只要求一阶导数存在 所以 FEM可以用 分段线形函数

弱形式把“点点都严格满足微分方程”变成“整体积分意义上满足方程”

Galerkin 方法：测试函数也选形函数

近似解：

$$ u_h(x)=\sum_j U_jN_j(x) $$ 弱形式：

$$ \int u_h'v'\,dx=\int fv\,dx $$ FEM/Galerkin 方法选择：

$$ v=N_i $$ 也就是说，用每个形函数当测试函数。

代入：

$$ \int \left( \sum_j U_jN_j' \right) N_i'\,dx = \int fN_i\,dx $$ 把求和拿出来：

$$ \sum_j U_j \int N_j'N_i'\,dx = \int fN_i\,dx $$ 定义矩阵：

$$ K_{ij} = \int N_i'N_j'\,dx $$ 定义右端：

$$ F_i = \int fN_i\,dx $$ 于是得到：

$$ \boxed{ K\mathbf U=\mathbf F } $$ 这就是 FEM 的矩阵方程。

K就是刚度矩阵 对于局部 就是 局部的刚度矩阵

组装就是：每个小单元算出一个局部矩阵，然后按节点编号加到全局矩阵的对应位置。

二维就用更高级的三角形

FEM 的基本步骤

1. 剖分区域：将求解域分割成大量小单元（二维用三角形，三维用四面体）
2. 选择插值函数：在每个单元内用多项式插值近似场量
3. 弱形式：将 PDE 乘以检验函数并分部积分，降低微分阶数
4. 组装与求解：叠加所有单元的贡献，得到线性方程组

例子

考虑：

$$ -\nabla^2\phi=f $$ 区域：

$$ \Omega $$ 边界：

$$ \partial\Omega $$ 假设 Dirichlet 边界：

$$ \phi=g \quad \text{on }\partial\Omega $$ 选一个测试函数 $v$，满足在 Dirichlet 边界上：

$$ v=0 $$ 将 PDE 乘以 $v$，积分：

$$ \int_\Omega (-\nabla^2\phi)v\,d\Omega = \int_\Omega fv\,d\Omega $$ 对左边分部积分。

利用恒等式：

$$ \nabla\cdot(v\nabla\phi) = \nabla v\cdot\nabla\phi + v\nabla^2\phi $$ 所以：

$$ v\nabla^2\phi = \nabla\cdot(v\nabla\phi) - \nabla v\cdot\nabla\phi $$ 因此：

$$ -\int_\Omega v\nabla^2\phi\,d\Omega = -\int_\Omega \nabla\cdot(v\nabla\phi)d\Omega + \int_\Omega \nabla v\cdot\nabla\phi\,d\Omega $$ 用散度定理：

$$ \int_\Omega\nabla\cdot(v\nabla\phi)d\Omega = \int_{\partial\Omega}v\frac{\partial\phi}{\partial n}dS $$ 所以：

$$ -\int_\Omega v\nabla^2\phi\,d\Omega = \int_\Omega\nabla v\cdot\nabla\phi\,d\Omega - \int_{\partial\Omega} v\frac{\partial\phi}{\partial n}dS $$ 对于 Dirichlet 边界，测试函数 $v=0$，边界项消失。

于是弱形式为：

$$ \boxed{ \int_\Omega\nabla v\cdot\nabla\phi\,d\Omega = \int_\Omega fv\,d\Omega } $$ 此为 FEM的基础 降低微分阶数 —— 因为有限元形函数通常是分片多项式，只要求一阶导存在即可！

Delaunay 三角剖分

有限元法中网格质量直接影响精度。Delaunay 三角剖分是最常用的剖分策略：

• ✅ 最大化所有三角形的最小角（避免"银针三角形"）
• ✅ 所有三角形的并集 = 点集的凸包
• ✅ 最近邻点必然由一条边相连

Poisson 方程 Dirichlet 问题的 FEM 推导

在三角形单元内，场量 $\phi$ 用三个顶点值插值：

$$ \phi(x,y) = N_a(x,y)\,\phi_a + N_b(x,y)\,\phi_b + N_c(x,y)\,\phi_c $$

其中 $N_a, N_b, N_c$ 是形函数（shape functions）。

通过变分法最小化能量泛函，得到局部刚度矩阵：

$$ k_{ab} = \int_\Delta \nabla N_a \cdot \nabla N_b \, dA $$

组装所有单元的局部刚度矩阵得到全局方程：

$$ K\phi = r $$

其中 $K$ 是对称矩阵，$r$ 是右端向量。

每个单元的局部刚度矩阵（具体形式）

设三角形三个顶点坐标为 $(x_a, y_a), (x_b, y_b), (x_c, y_c)$：

def local_stiffness(verts):
    """计算单个三角形单元的局部刚度矩阵"""
    xa, ya = verts[0]
    xb, yb = verts[1]
    xc, yc = verts[2]
    # 面积的两倍
    J = (xb - xa)*(yc - ya) - (xc - xa)*(yb - ya)
    # 形函数梯度分量
    bx = (yb - yc) / J  # ∂Na/∂x
    cx = (xc - xb) / J  # ∂Na/∂y → 对应 Na
    # ... 轮换得到其他顶点的分量
    # 局部刚度矩阵 3×3
    k = np.zeros((3, 3))
    k[0, 0] = (bx*bx + cx*cx) * J / 2  # 依此类推
    return k

回顾展开推导

三角形线性单元

设一个三角形单元 $\Delta$ 有三个顶点：

$$ a,b,c $$ 坐标分别为：

$$ (x_a,y_a),\quad (x_b,y_b),\quad (x_c,y_c) $$ 在这个三角形内部，用线性插值：

$$ \phi(x,y) = N_a(x,y)\phi_a + N_b(x,y)\phi_b + N_c(x,y)\phi_c $$ 其中 $N_a,N_b,N_c$ 是形函数。

它们满足：

$$ N_a(a)=1,\quad N_a(b)=0,\quad N_a(c)=0 $$ 类似地：

$$ N_b(b)=1,\quad N_c(c)=1 $$ 并且：

$$ N_a+N_b+N_c=1 $$

线性三角形形函数具体形式

形函数可以写成：

$$ N_i(x,y) = \frac{a_i+b_ix+c_iy}{2A} $$ 其中 $A$ 是三角形面积。

对顶点 $a$：

$$ a_a=x_by_c-x_cy_b $$ 对顶点 $b$：

$$ a_b=x_cy_a-x_ay_c $$ 对顶点 $c$：

$$ a_c=x_ay_b-x_by_a $$ 面积满足：

$$ 2A = \begin{vmatrix} 1&x_a&y_a\\ 1&x_b&y_b\\ 1&x_c&y_c \end{vmatrix} $$ 也就是：

$$ 2A = (x_b-x_a)(y_c-y_a) - (x_c-x_a)(y_b-y_a) $$ 通常取绝对值保证面积为正。

形函数梯度

因为 $N_i$ 是线性的，所以梯度是常数：

$$ \nabla N_i = \left( \frac{\partial N_i}{\partial x}, \frac{\partial N_i}{\partial y} \right) = \left( \frac{b_i}{2A}, \frac{c_i}{2A} \right) $$

局部刚度矩阵推导

弱形式对应的局部刚度矩阵为：

$$ k_{ij}^{(e)} = \int_{\Delta_e} \nabla N_i\cdot\nabla N_j\,dA $$ 由于 $\nabla N_i$ 和 $\nabla N_j$ 在三角形内是常数：

$$ k_{ij}^{(e)} = \left( \nabla N_i\cdot\nabla N_j \right) \int_{\Delta_e}dA $$ 所以：

$$ \boxed{ k_{ij}^{(e)} = \frac{b_ib_j+c_ic_j}{4A} } $$ 这就是线性三角形单元的局部刚度矩阵。

局部载荷向量

右端：

$$ r_i^{(e)} = \int_{\Delta_e} f N_i\,dA $$ 如果 $f$ 在单元内近似为常数，可以取三角形重心处值：

$$ f_e=f(x_{\text{centroid}},y_{\text{centroid}}) $$ 则：

$$ r_i^{(e)} \approx f_e\int_{\Delta_e}N_i\,dA $$ 线性三角形中：

$$ \int_{\Delta_e}N_i\,dA=\frac{A}{3} $$ 所以：

$$ \boxed{ r_i^{(e)} \approx \frac{A}{3}f_e } $$

组装全局矩阵

每个单元有一个 $3\times 3$ 局部刚度矩阵：

$$ k^{(e)} $$ 假设该单元三个局部节点对应全局编号：

$$ I_a,I_b,I_c $$ 那么把局部矩阵加到全局矩阵：

$$ K_{I_p I_q} \mathrel{+}= k_{pq}^{(e)} $$ 右端同理：

$$ r_{I_p}\mathrel{+}=r_p^{(e)} $$ 最后得到：

$$ \boxed{ K\boldsymbol\phi=\mathbf r } $$

Dirichlet 边界怎么处理？

如果某个边界点已知：

$$ \phi_i=g_i $$ 常见处理方式：

1. 把该自由度从未知量中删除；
2. 或者修改矩阵第 $i$ 行：

$$ K_{ii}=1,\quad K_{ij}=0\ (j\neq i) $$ 右端设为：

$$ r_i=g_i $$ 为了保持对称性，通常还要同时修改列。

Neumann 边界在 FEM 中很自然

如果边界给：

$$ \frac{\partial\phi}{\partial n}=q $$ 弱形式中边界项不会消失：

$$ -\int_{\partial\Omega} v\frac{\partial\phi}{\partial n}dS $$ 根据方程符号约定，它会进入右端。

所以 FEM 里 Neumann 条件是自然边界条件。

FEM 本征值问题中的质量矩阵

如果求 Schrödinger 方程：

$$ H\psi=E\psi $$ FEM 展开：

$$ \psi=\sum_i c_iN_i $$ 投影后得到：

$$ H\mathbf c=E M\mathbf c $$ 这里 $M$ 是质量矩阵：

$$ M_{ij} = \int N_iN_j\,d\Omega $$ 对线性三角形，局部质量矩阵为：

$$ M^{(e)} = \frac{A}{12} \begin{pmatrix} 2&1&1\\ 1&2&1\\ 1&1&2 \end{pmatrix} $$ 这和 B-样条里的重叠矩阵 $S$ 是同一个角色。

FEM vs 有限差分法

第十章：边界条件进阶 — 多极展开、周期性、Twisted

💡 物理直觉

前面讨论的是"给定边界值"的简单情况。实际科研中，边界条件本身就是个难题。比如 Coulomb 势在无穷远处衰减极慢（~1/r），怎么截断？

1. Coulomb 势边界条件 — 多极展开法

对于长程 Coulomb 势，在边界外面还有贡献。处理方法：

步骤：将边界上的源 $\rho_b$ 分离出来：

$$ \nabla^2\phi = -4\pi\rho $$

在边界上的贡献用多极展开（Spherical Harmonics $Y_{\ell m}$）解析计算：

$$ \phi_{\text{boundary}} = \sum_{\ell,m} \frac{Q_{\ell m}}{r^{\ell+1}} Y_{\ell m}(\theta, \phi) $$

其中 $Q_{\ell m} = \int r^\ell Y_{\ell m}^*(\theta, \phi) \,\rho(\mathbf{r})\, d^3r$ 是多极矩。

然后用有限差分法求解内部：$\nabla^2\phi = -4\pi(\rho - \rho_b)$（减去边界源的贡献）。

说人话：把长程部分用球谐函数解析算了，剩下的短程部分用差分法。分而治之！

2. 周期性边界条件 PBC

傅里叶微分算符天然处理周期性边界：$\phi(x+L) = \phi(x)$。

在固体物理、晶格 QCD 中，PBC 是标配。

但有一个问题：有限盒子的假量子化效应——本该连续的谱被离散化了。

3. Twisted 边界条件 TBC

Twisted 边界条件是对 PBC 的推广：

$$ \phi(x+L) = e^{i\theta} \phi(x) $$

• 物理动机：通过变化 twist 角 $\theta$ 来恢复连续谱
• 实现：修改傅里叶动量 $k = (2\pi n + \theta)/L$
• 应用：晶格 QCD、量子 Monte Carlo、凝聚态物理、中子星壳层 Pasta 结构

4. 吸收边界条件 ABC

对于散射态/连续谱问题，不让波函数在边界反射回来：

入射波 → 计算区域 → 透射（吸收）
           ↛ 反射回来！

三种边界条件对比

第十一章：Jacobi 松弛迭代法

💡 物理直觉

对于 Poisson 方程 $\nabla^2 u = -f$，离散化后得到 $A\mathbf{u} = \mathbf{b}$。直接求逆对大规模矩阵不现实（O(N³)）。Jacobi 迭代的一条思路："方程两边同时加回主对角项，变成不动点迭代。"

先乱猜一个解， 然后反复使用“邻居平均值”进行修正，直到整体平均平衡！！！

Jacobi 迭代公式推导

从五点差分公式出发：

$$ \frac{u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4u_{i,j}}{h^2} = -f_{i,j} $$

移项，解出 $u_{i,j}$：

$$ u_{i,j}^{\text{(new)}} = \frac{1}{4} \left[ u_{i+1,j}^{\text{(old)}} + u_{i-1,j}^{\text{(old)}} + u_{i,j+1}^{\text{(old)}} + u_{i,j-1}^{\text{(old)}} + h^2 f_{i,j} \right] $$

说人话：每个格点的新值 = 四个邻居旧值的平均值 + 源项修正。不断重复直到收敛。本质上，这就是让势函数慢慢"松弛"到平衡态！它会慢慢把不平衡的势函数抹平，直到满足平衡条件。

数值实现

import numpy as np

def jacobi_poisson_2d(u, f, h, tol=1e-6, max_iter=10000):
    u_old = u.copy()
    u_new = u.copy()

    for it in range(max_iter):
        u_old[:, :] = u_new

        u_new[1:-1, 1:-1] = 0.25 * (
            u_old[2:, 1:-1]
            + u_old[:-2, 1:-1]
            + u_old[1:-1, 2:]
            + u_old[1:-1, :-2]
            + h**2 * f[1:-1, 1:-1]
        )

        diff = np.max(np.abs(u_new - u_old))

        if diff < tol:
            print("Jacobi 收敛步数:", it)
            break

    return u_new

加速技巧

• Gauss-Seidel 迭代：使用最新值而非旧值（Jacobian→Gauss-Seidel，收敛翻倍）
• 超松弛 SOR：加一个加速因子 $\omega$（$1 < \omega < 2$）
• 多重网格法：粗网格上快速收敛低频分量，细网格上打磨高频分量

总结与对比

核心思想归纳

1. 多维 PDE 的密钥 → 投影转换为一维：矩阵形式的本质就是把 2D/3D 格点压扁成一个大向量
2. 不同算法的精度关键 → 微分算符：如何逼近导数决定了方法的上限
3. 两大流派：基空间（灵活高效）vs 格点空间（精确但计算量大）

七种方法终极对比表

📊 方法选择决策树

遇到 PDE 问题
├── 规则区域 + 周期边界 → FFT 谱方法（最快）
├── 规则区域 + 一般边界 → 有限差分法（最直接）
├── 束缚态本征值问题
│   ├── 基函数合适 → 基矢展开法（最高效）
│   ├── 需要高精度 → DVR（最优收敛）
│   └── 需要非均匀网格 → B-样条（最灵活）
├── 多尺度 / 自适应需求 → 多小波分析（最智能）
├── 不规则复杂边界 → 有限元法（最通用）
└── 超大规模稀疏 → Jacobi / 多重网格迭代

🔗 PDE 方法与物理问题的对应

💡 课程金句备忘

• "多维问题的关键是投影转换为一维"
• "不同算法的精度关键是微分算符"
• "基空间灵活高效，格点空间精确但计算量大"
• "利用对称性和好量子数可以降维，打破对称性则计算量增加几个量级"
• "Nyquist 准则：格点间距不能大于最小波长的 1/2，否则信息丢失"

笔记结束 · 2026-06-07