求解泊松方程：有限差分法在梯形域的二维网格映射方法咨询

嘿，很高兴你已经搞定了规则网格上的泊松方程有限差分实现！梯形域其实是结构化非规则域里比较好处理的一类，核心思路就是你提到的坐标映射——把梯形域（我们叫它「物理域$\Omega$」）和一个单位正方形/规则矩形（「计算域$\hat{\Omega}$」，比如$[0,1]\times[0,1]$）建立一一对应的双向映射，然后在计算域上用你熟悉的规则网格做有限差分，最后把结果映射回物理域就行。

下面一步步给你拆解具体怎么做：

一、构建梯形域与规则计算域的双向映射 #

先假设你的梯形是最常见的那种：下底在$y=0$，从$x=a$到$x=b$；上底在$y=H$，从$x=c$到$x=d$（$a < c < d < b$或者反过来都没关系，调整参数就行）。

选项1：仿射线性映射（最简单，优先用） #

这种映射适合所有边都是直线的梯形，完全线性，计算成本最低，推导也简单。我们定义物理坐标$(x,y)$和计算域坐标$(\xi,\eta)$（$\xi\in[0,1], \eta\in[0,1]$）的关系：

• $y$方向直接对应：$\eta$从0到1，对应$y$从0到$H$，所以 $y = H\eta$
• $x$方向是随$\eta$线性变化的：每个高度$\eta$对应的$x$范围，左边界从$a$（$\eta=0$）线性变到$c$（$\eta=1$），右边界从$b$线性变到$d$。所以：
  $x_L(\eta) = a + (c - a)\eta$（左边界的$x$随$\eta$的变化）
  $x_R(\eta) = b + (d - b)\eta$（右边界的$x$随$\eta$的变化）
  然后任意$\xi$对应的$x$就是：$x = x_L(\eta) + (x_R(\eta) - x_L(\eta))\xi$

反过来，从物理域转计算域的逆映射也很容易推：

• $\eta = y/H$
• $\xi = \frac{x - x_L(\eta)}{x_R(\eta) - x_L(\eta)}$

这种映射的好处是：雅可比矩阵的元素要么是常数，要么只和$\eta$有关，后续推导泊松方程会省很多事。

选项2：双线性参数化映射（更灵活） #

如果你的梯形有轻微的曲线边（或者是近似梯形的曲线域），可以用双线性映射，公式是：

x = Aξ + Bη + Cξη + D
y = Eξ + Fη + Gξη + H

把梯形的四个顶点坐标代入，解出8个系数A-H就行。比如四个顶点是$(a,0), (b,0), (d,H), (c,H)$，代入后解线性方程组就能得到所有系数，这种映射也是双向可逆的，不过对于标准直线梯形，仿射映射足够了。

二、推导变换后的泊松方程 #

原物理域的泊松方程是：
$$\nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f(x,y)$$
其中$u$是待求函数，$f$是源项。

我们需要把这个方程转换到计算域$(\xi,\eta)$下，核心是用链式法则求偏导。先写一阶偏导的关系：
$$\frac{\partial u}{\partial \xi} = \frac{\partial u}{\partial x}\frac{\partial x}{\partial \xi} + \frac{\partial u}{\partial y}\frac{\partial y}{\partial \xi}$$
$$\frac{\partial u}{\partial \eta} = \frac{\partial u}{\partial x}\frac{\partial x}{\partial \eta} + \frac{\partial u}{\partial y}\frac{\partial y}{\partial \eta}$$

把这个写成矩阵形式，就是：
$$\begin{pmatrix} \partial u/\partial \xi \ \partial u/\partial \eta \end{pmatrix} = J \begin{pmatrix} \partial u/\partial x \ \partial u/\partial y \end{pmatrix}$$
这里的$J$就是雅可比矩阵：
$$J = \begin{pmatrix} \partial x/\partial \xi & \partial y/\partial \xi \ \partial x/\partial \eta & \partial y/\partial \eta \end{pmatrix}$$

再推导二阶偏导，最终变换到计算域的泊松方程会变成：
$$\frac{1}{\det(J)} \left[ \frac{\partial}{\partial \xi} \left( \frac{J_{22}\frac{\partial u}{\partial \xi} - J_{12}\frac{\partial u}{\partial \eta}}{\det(J)} \right) + \frac{\partial}{\partial \eta} \left( \frac{-J_{21}\frac{\partial u}{\partial \xi} + J_{11}\frac{\partial u}{\partial \eta}}{\det(J)} \right) \right] = f(x(\xi,\eta), y(\xi,\eta))$$
其中$\det(J) = (\partial x/\partial \xi)(\partial y/\partial \eta) - (\partial x/\partial \eta)(\partial y/\partial \xi)$是雅可比行列式，$J_{ij}$是雅可比矩阵的第i行第j列元素。

重点：仿射映射下的简化 #

对于我们之前的仿射线性映射，雅可比矩阵的元素计算很简单：

• $\partial x/\partial \xi = x_R(\eta) - x_L(\eta)$（只和$\eta$有关）
• $\partial x/\partial \eta = (c - a) + (d - b - c + a)\xi$
• $\partial y/\partial \xi = 0$
• $\partial y/\partial \eta = H$

所以雅可比行列式$\det(J) = H \times (x_R(\eta) - x_L(\eta))$，也是只和$\eta$有关的函数，代入变换后的方程会简化很多，计算量大大降低。

三、在计算域上用有限差分求解 #

现在计算域是规则网格，比如我们取$N\times M$的网格，$\xi_i = i/(N-1)$（$i=0,1,...,N-1$），$\eta_j = j/(M-1)$（$j=0,1,...,M-1$）。

步骤如下：

1. 预计算网格点映射：对每个计算域网格点$(\xi_i,\eta_j)$，用之前的映射公式算出对应的物理域坐标$(x_{ij}, y_{ij})$，同时预计算每个点的雅可比矩阵元素和行列式。
2. 离散变换后的方程：用你熟悉的有限差分格式（比如中心差分）近似偏导数。比如二阶偏导的中心差分：
   $$\frac{\partial^2 u}{\partial \xi^2}|{i,j} \approx \frac{u{i+1,j} - 2u_{i,j} + u_{i-1,j}}{(\Delta\xi)^2}$$
   同理处理$\eta$方向的偏导，把变换后的方程里的每一项都用差分近似替换。
3. 组装线性方程组：把离散后的方程整理成$Au = b$的形式，其中$A$是系数矩阵，$u$是未知量向量，$b$是右端项（包含源项$f$和边界条件）。
4. 求解与映射回物理域：用你之前会的求解器（比如Gauss-Seidel、共轭梯度法）解线性方程组，得到计算域上的$u_{ij}$，再用逆映射把每个$u_{ij}$对应到物理域的$(x_{ij}, y_{ij})$上，就是最终结果。

四、边界条件的处理 #

边界条件也要同步转换：

• Dirichlet边界：如果物理域上是$u(x,y)=g(x,y)$，那么计算域的边界点（$\xi=0, \xi=1, \eta=0, \eta=1$）对应的$u$值就是$g(x(\xi,\eta), y(\xi,\eta))$，直接代入就行。
• Neumann边界：如果是$\frac{\partial u}{\partial n} = h(x,y)$，需要把物理域的法向导数转换到计算域上，用链式法则推导：先算出物理域边界的法向量，再通过雅可比矩阵的逆转换为计算域的向量，然后计算对应的导数近似。

小技巧：验证映射与求解的正确性 #

• 先验证映射：取梯形的四个顶点，用映射公式计算，看是否对应计算域的四个顶点$(0,0),(1,0),(1,1),(0,1)$，确保双向映射是正确的。
• 用解析解测试：找一个有解析解的泊松方程（比如$u(x,y)=x^2 + y^{2$，此时$\nabla}2 u=4$，源项$f=4$），用你的方法求解，对比计算结果和解析解，验证整个流程的正确性。

内容的提问来源于stack exchange，提问作者Helena Coquand