核心主题：有限差分法、基矢展开法（含 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) 通过函数耦合在一起。

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

从一维到二维的推广

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

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

五点查分的具体推导：

二维 Laplace 算子是：

在均匀网格上：

函数值记为：

对 $x$ 方向二阶导：

对 $y$ 方向二阶导：

两者相加：

这就是二维五点差分格式。

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

若：

则：

整理：

如果 $h_x=h_y=h$，就退化为五点公式。

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

五点公式：

可以写成：

也就是说 约等于 邻居的平均值 减去 自身 就是 $\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倍自己，再除以格距平方。这本质上是在问："我和邻居们的平均值差了多少？"

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

即，没有源的时候，一个点的值等于周围“邻居” 的平均值

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

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

💡 物理直觉

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

重要提醒：两种矩阵约定

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

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

写成：

其中 $A_0$ 的主对角是 $4$，邻居位置是 $-1$。

此时右端应该是：

约定 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}$。按行（或列）排列：

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

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

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

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

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

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

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

离散方程：

其中：

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

把边界已知量移到右端：

所以边界值是以正号加到 RHS 上的。

一般规律：

这是在矩阵不带 $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矩阵！

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

一维二阶差分矩阵： $$ TN= \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( Iy\otimes Tx + Ty\otimes Ix \right) $$ 如果 $Nx=N_y=3$，就会得到你写的 $9\times 9$ 块三对角矩阵。

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

矩阵性质

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

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

对称：

正定：

只要 $\mathbf u\neq 0$。

正定性来自物理能量：

离散后对应：

所以这个矩阵适合用共轭梯度法 CG

📊 矩阵规模与计算代价

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

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

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

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

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

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

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

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

否则解通常不唯一

三类边界条件对比

Neumann 边界条件的数值处理

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

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

一维 Neumann 的 虚格点 推导

考虑一维 Poisson：

左边界给 Neumann：

在 $x=0$ 处引入虚点：

中心差分近似一阶导数：

令它等于 $g$：

解出虚点：

现在在边界点 $x_0$ 上离散 PDE：

二阶导：

所以：

代入：

得到：

整理：

因此：

这就是 Neumann 边界点对应的离散方程。

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

则：

这意味着边界外的虚点是边界内点的镜像。

二维 Neumann 边界

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

引入虚点 $u_{-1,j}$，有：

所以：

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

写成：

代入虚点：

得到：

这就把虚点消掉了。

Robin 边界条件

Robin 条件：

例如一维左边界：

用中心差分：

代入：

解出虚点：

即：

然后再代回 PDE 的差分方程。

纯Neumann问题的特殊性

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

则解只确定到一个常数。

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

也是解。

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

还需要兼容条件。

由散度定理：

因为：

所以：

如果这个条件不满足，纯 Neumann 问题无解。

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

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

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

💡 物理直觉

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

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

在格点空间做微分的"近邻相减"，在动量空间就是简单的乘法！

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

离散 Fourier变换 DFT

在周期区间：

离散 Fourier 展开：

其中波数：

但是 DFT 的频率顺序通常是：

在 NumPy 中：

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

会自动给出正确频率。

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

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

Fourier 微分矩阵

令 DFT 矩阵为 $F$，则：

反变换：

一阶导数在 Fourier 空间是乘 $ik$，所以：

回到实空间：

因此一阶微分矩阵是：

二阶微分矩阵：

因为：

所以：

注意：这个关系在 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 微分

二维周期函数：

则：

在 Fourier 空间变成：

所以算法是：

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方法 默认周期边界条件！！！

也就是说：

如果函数不周期，强行用 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 方程

考虑周期边界下：

Fourier 空间：

所以：

但当：

分母为零。

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

必须要求：

也就是：

然后可以人为设定：

即选择零平均势。

FFT 与共振展宽

如果一个信号：

它的 Fourier 变换近似为：

整理指数：

积分：

谱强度：

这就是 Lorentz 线形。

越大，谱线越宽，寿命越短：

若带 $\hbar$，则：

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

💡 物理直觉

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

基本公式推导

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

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

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

写成矩阵形式：

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

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

随着基函数数目增大：

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

3D 笛卡尔谐振子基

三维笛卡尔谐振子基：

其中一维谐振子基函数：

$H_n$ 是 Hermite 多项式，$b$ 是振子长度参数。

总壳层量子数：

如果截断到：

则基函数数目是：

这个数目随着 $N_{\max}$ 迅速增长。

轴对称体系（柱坐标）

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

基函数可以分离为：

其中：

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

常用截断：

保留：

好量子数为什么能降维？

如果 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$，格距：

对于无限 sinc-DVR，常见二阶导数矩阵是：

于是动能：

对应：

周期 sinc-DVR 或 Fourier-DVR 的公式会略有不同，比如会出现：

这正是你代码里的结构。

需要注意：

不同文献里 $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

给定正交基函数：

定义投影算符：

这个投影算符的意思是： 我们只在这 $N$ 个基函数张成的有限维空间里面近似真实波函数。

选择高斯积分点和权重

选一组积分点：

以及对应权重：

满足高斯积分近似：

你讲义里写的是 gauge points，更准确地说一般叫：

或者：

如果这组点和权重选得好，那么对于基函数乘积：

可以用求和精确或高精度表示：

DVR 变换矩阵

定义矩阵：

这里：

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

因为：

所以：

也就是：

在理想情况下，$U$ 是酉矩阵或正交矩阵。

DVR 基函数的构造

DVR 基函数定义为：

如果基函数都是实函数，可以省略复共轭：

写成坐标形式：

这个 $\chi_\alpha(x)$ 就是 DVR 基函数。

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

我们希望它满足类似：

也就是：

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

来推导。

在格点 $x_\beta$ 处：

利用矩阵 $U$ 的正交性：

代入：

所以：

两边除以 $\sqrt{w_\beta}$：

左边正是 $\chi_\alpha(x_\beta)$。

因此：

如果重新定义归一化后的 DVR cardinal 函数：

那么：

当 $\alpha=\beta$ 时：

当 $\alpha\neq\beta$ 时：

所以：

这就是 DVR 的核心插值性质。

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

考虑势能矩阵：

写成积分：

用高斯积分近似：

利用：

以及：

代入：

整理：

只有当：

时有贡献。

所以：

这就是 DVR 最重要的结果。

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

动能矩阵如何处理？

动能算符为：

在 FBR 中，动能矩阵是：

这通常可以解析算出，或者用高精度积分算出。

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

即：

所以 Hamiltonian 矩阵变成：

最终求解：

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-样条的阶数；

多项式次数为：

例如：

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

一阶 B-样条定义为：

高阶通过递推构造：

如果分母为零，则对应项定义为零。

由递推关系可以看出：

只在区间：

上非零。

所以它是紧支撑的。

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

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

因此重叠矩阵和 Hamiltonian 矩阵都是带状矩阵。

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

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

将波函数展开：

代入方程：

利用线性性：

左乘测试函数 $B_i^M(x)$，并积分：

交换求和与积分：

定义：