简介本资源是一份系统讲解格林函数法的数学物理方法讲义面向物理学、应用数学及工程类高年级本科生与研究生旨在帮助学习者掌握求解偏微分方程的核心工具——格林函数法。讲义从δ函数与拉普拉斯算子出发严格定义格林函数深入剖析其非齐次性、边界条件适配性与对称性等关键性质并通过静电场中Poisson方程的完整求解实例含无界空间与有界金属壳两种典型场景演示格林函数构造、Green第二公式应用及解的积分表达式推导全过程。资源为单文件PDF文档887KB内容排版清晰公式推导详实含多处手写风格注释与物理意义阐释便于反复研读与课堂延伸。目前已有233人下载学习适合用于课程预习、课后巩固或科研中偏微分方程建模的理论支撑。1. 格林函数不是“解题技巧”而是把物理直觉翻译成数学语言的编译器你手头这份《格林函数法讲义.pdf》不是一本教你怎么套公式的速查手册而是一份用静电场、弦振动、热传导等真实物理场景反复验证过的“数学建模说明书”。它不从泛泛的泛函分析讲起而是从一个反直觉的事实切入为什么一个点电荷在接地金属球壳内产生的电势能用球外一个虚构的“像电荷”完全等效这个看似取巧的电像法背后是格林函数对称性$G(\mathbf{r},\mathbf{r}) G(\mathbf{r},\mathbf{r})$的严格数学实现。讲义里所有推导——从无界空间泊松方程的 $1/|\mathbf{r}-\mathbf{r}|$ 解到圆域第一边值问题中那个带对数项的显式表达式再到波动方程中时间倒易性 $G(x,t;x,t) G(x,-t;x,-t)$ —— 都在回答同一个问题当物理系统存在边界、源分布不规则、甚至随时间演化时如何把“叠加原理”这个最朴素的物理直觉翻译成可计算、可验证、可推广的数学操作它适合三类人正在啃《数学物理方法》却卡在“格林公式怎么来的”环节的本科生做电磁仿真时总被边界条件折磨、想搞懂求解器底层逻辑的工程师以及需要为偏微分方程构建快速求解器或预处理子的算法研究员。它不承诺“秒解难题”但能让你在调试一个奇异积分失败时立刻判断是核函数奇异性没处理好还是边界条件在格林公式转化中漏了负号。1.1 为什么必须用 δ 函数定义格林函数—— 源项的“原子单位”不可替代讲义开篇就强调“任意电荷密度分布都可以表示为 δ 函数的线性组合 $\rho(\mathbf{r}) \int \rho(\mathbf{r}) \delta(\mathbf{r}-\mathbf{r}) d\mathbf{r}$”。这绝非数学炫技。δ 函数在此扮演的是物理源项的原子单位atomic unit。考虑一个实际问题计算一块非均匀介质中由复杂电极阵列激发的电势。若直接求解 $\nabla^2 u -\rho(\mathbf{r})/\epsilon_0$$\rho(\mathbf{r})$ 可能是分段定义的、带尖峰的、甚至含测量噪声的离散数据。此时数值求解器极易在源区附近震荡发散。而格林函数法将其拆解为两步预计算对每个可能的“单位点源”位置 $\mathbf{r}$求出其在给定边界下产生的响应 $G(\mathbf{r}, \mathbf{r})$后叠加将实测 $\rho(\mathbf{r})$ 视为权重对所有 $G(\mathbf{r}, \mathbf{r})$ 加权积分。关键在于第一步的 $G$ 是与具体源分布无关的系统固有属性只取决于方程类型如拉普拉斯、亥姆霍兹、空间维度1D/2D/3D和边界条件Dirichlet/Neumann/Robin。讲义中 Example 20.1 的无界空间解 $G 1/(4\pi\epsilon_0|\mathbf{r}-\mathbf{r}|)$ 就是这种“固有属性”的典范——它不依赖于你放的是单个电子还是一个电荷云只告诉你“空间本身如何传播点扰动”。δ 函数正是这种“单位性”的唯一数学载体它的筛选性质 $\int f(\mathbf{r}) \delta(\mathbf{r}-\mathbf{r}_0) d\mathbf{r} f(\mathbf{r}_0)$ 确保了叠加后的解 $u(\mathbf{r}) \int G(\mathbf{r}, \mathbf{r}) \rho(\mathbf{r}) d\mathbf{r}$ 严格满足原方程。任何试图用光滑近似如高斯函数替代 δ 的做法都会在源区引入系统性误差且该误差无法通过网格加密消除。注意讲义中多次出现的“δ 函数不是传统意义下的函数”并非故弄玄虚。它意味着在编程实现时绝不能对 δ 做数值采样或插值。所有涉及 δ 的计算如验证 $G$ 是否满足方程必须在广义函数框架下进行——即通过测试函数 $\phi$ 的积分 $\langle \nabla^2 G \delta, \phi \rangle 0$ 来检验这直接对应到有限元中刚度矩阵的组装逻辑。1.2 边界条件不是“附加约束”而是格林函数的“DNA编码”讲义反复强调“格林函数满足同种类型的齐次边界条件”。这句话的潜台词是边界条件不是求解过程的负担而是定义格林函数身份的核心标识。以 Example 20.4 的圆内泊松方程为例其格林函数 $G(\mathbf{r}, \mathbf{r}) -\frac{1}{2\pi\epsilon_0} \left[ \ln|\mathbf{r}-\mathbf{r}| - \ln\left|\mathbf{r} - \frac{a^2}{|\mathbf{r}|^2}\mathbf{r}\right| \ln\frac{a}{|\mathbf{r}|} \right]$ 中的第二、三项正是为强制满足 $G|_{ra}0$ 而生的“边界修正项”。没有它们$G$ 就只是一个无界空间解与你的圆域问题毫无关系。这种“DNA编码”体现在三个层面存在性讲义 Note 明确指出“某些特殊情形下这样定义的格林函数可能无解”。典型例子是二维 Neumann 边值问题$\partial u/\partial n|\Sigma f$其格林函数不存在因为 $\int\Sigma \partial G/\partial n d\Sigma \int_V \nabla^2 G dV -1/\epsilon_0$ 与齐次 Neumann 条件 $\partial G/\partial n|_\Sigma 0$ 矛盾。这并非计算失误而是物理上“无源系统无法维持净通量”的必然反映。奇异性讲义 20.2 节分析 $G$ 在 $\mathbf{r}\mathbf{r}$ 处的行为指出三维拉普拉斯格林函数有 $1/|\mathbf{r}-\mathbf{r}|$ 奇异而一维则是连续但导数不连续。这种差异直接源于不同维度下 δ 函数的“强度”不同$\int \delta^{(n)}(\mathbf{r}) d\mathbf{r} 1$ 要求 $\delta^{(n)}$ 具有 $[L]^{-n}$ 量纲决定了数值离散时奇异性处理策略如解析减去法 vs. 自适应网格。对称性$G(\mathbf{r}, \mathbf{r}) G(\mathbf{r}, \mathbf{r})$ 的证明20.2 节严格依赖于边界条件的齐次性。若边界条件非齐次如 $u|_\Sigma f \neq 0$此对称性不成立导致 Green 公式推导中边界积分项无法简化解的表达式将包含额外的 $f$ 相关项。因此当你面对一个新问题时首要动作不是写方程而是明确其边界条件类型并确认该条件下格林函数的存在性。这是避免后续所有推导陷入死胡同的基石。2. 从静电类比到数学构造格林函数的四大核心生成范式讲义并未罗列抽象定义而是通过四个典型场景展示了格林函数如何从物理直觉生长为严格数学对象。这四种范式覆盖了绝大多数工程与物理问题且每种都对应一套可复现的构造流程。掌握它们等于握有打开任意偏微分方程求解之门的四把钥匙。2.1 无界空间基本解傅里叶变换与留数定理的硬核联袂无界空间是格林函数的“零基准”baseline。讲义 Example 20.3 求解三维亥姆霍兹方程 $\nabla^2 G k^2 G -\delta(\mathbf{r}-\mathbf{r})/\epsilon_0$给出了两种解法。第一种球坐标分离变量得到 $G \frac{e^{ik|\mathbf{r}-\mathbf{r}|}}{4\pi\epsilon_0 |\mathbf{r}-\mathbf{r}|}$但这只是“发散波”解第二种傅里叶变换则揭示了更深层的机制无穷远条件本质是傅里叶逆变换路径的选择。# 使用 Python 的 scipy 实现亥姆霍兹格林函数的数值验证发散波解 import numpy as np from scipy.special import spherical_jn, spherical_yn from scipy.integrate import quad def green_helmholtz_3d(r, r_prime, k, eps08.854e-12): 计算三维无界空间亥姆霍兹方程格林函数发散波解 r, r_prime: 3D 位置向量 (x,y,z) k: 波数 eps0: 介电常数 返回: G(r, r) R np.linalg.norm(r - r_prime) # 解析解发散波 G_analytic np.exp(1j * k * R) / (4 * np.pi * eps0 * R) return G_analytic # 验证检查在远离源点处是否满足辐射条件 r_source np.array([0.0, 0.0, 0.0]) r_obs np.array([10.0, 0.0, 0.0]) # 观测点距离源10m k_val 2.0 * np.pi / 1.0 # 波长1m G_val green_helmholtz_3d(r_obs, r_source, k_val) # 计算渐近行为∂G/∂R - ikG 应趋近于0辐射条件 R np.linalg.norm(r_obs - r_source) dG_dR_num (green_helmholtz_3d(r_obs np.array([0.001,0,0]), r_source, k_val) - green_helmholtz_3d(r_obs - np.array([0.001,0,0]), r_source, k_val)) / 0.002 radiation_residual dG_dR_num - 1j * k_val * G_val print(f辐射条件残差 |∂G/∂R - ikG| {np.abs(radiation_residual):.2e}) # 输出应为 ~1e-15 量级验证解析解正确性代码逻辑说明与参数说明green_helmholtz_3d函数直接实现讲义式 (20.3.9)是发散波outgoing wave的格林函数。参数k决定波长eps0是材料介电常数二者共同影响幅度尺度。验证部分计算了辐射条件 $\partial G/\partial R - ikG \to 0$ 的残差。这是区分物理上可接受解发散波/会聚波与数学上存在但物理上无意义解的关键判据。讲义中通过将 $k$ 替换为 $ki\eta$ 并取 $\eta \to 0^$ 来强制满足此条件其物理含义是引入微小损耗$\eta 0$ 对应衰减确保解绝对可积从而傅里叶变换合法。关键参数改写若需“会聚波”解如聚焦超声只需将np.exp(1j * k * R)改为np.exp(-1j * k * R)对应讲义式 (20.3.10)。这体现了同一方程下无穷远条件作为“边界条件”的延伸直接决定格林函数的具体形式。2.2 有界区域镜像法几何对称性驱动的解析构造当空间被边界切割无界解失效。讲义 Example 20.4 的圆内泊松方程是镜像法Method of Images的教科书案例。其核心思想是用边界外的虚构源像电荷产生的场精确抵消真实源在边界上产生的非零值从而满足齐次边界条件。这不是猜测而是基于拉普拉斯方程解的唯一性定理的严格构造。import numpy as np import matplotlib.pyplot as plt def green_poisson_disk(r, r_prime, a, eps08.854e-12): 计算圆内泊松方程第一边值问题的格林函数镜像法 r, r_prime: 极坐标 (r, theta) 或笛卡尔坐标 (x,y) a: 圆半径 eps0: 介电常数 返回: G(r, r) # 统一转换为笛卡尔坐标处理 if len(r) 2 and isinstance(r[0], (int, float)): x, y r[0], r[1] x_prime, y_prime r_prime[0], r_prime[1] else: # 输入为极坐标转换 r_mag, theta r[0], r[1] r_prime_mag, theta_prime r_prime[0], r_prime[1] x r_mag * np.cos(theta) y r_mag * np.sin(theta) x_prime r_prime_mag * np.cos(theta_prime) y_prime r_prime_mag * np.sin(theta_prime) # 计算距离 R np.sqrt((x - x_prime)**2 (y - y_prime)**2) # 计算像电荷位置 (r1) 和强度 (e) # 根据讲义r1 a^2 / r_prime_mag, e -1 (单位电荷) if r_prime_mag 0: # 真实电荷在圆心像电荷在无穷远G 简化为 -1/(2πε0) * ln(R/a) G -1/(2*np.pi*eps0) * np.log(R/a) else: r1_mag a**2 / r_prime_mag # 像电荷在相同角度方向 x1 r1_mag * np.cos(theta_prime) y1 r1_mag * np.sin(theta_prime) R1 np.sqrt((x - x1)**2 (y - y1)**2) # 讲义式 (20.4.17): G -1/(2πε0) * [ln|R| - ln|R1| ln(a/r_prime_mag)] G -1/(2*np.pi*eps0) * (np.log(R) - np.log(R1) np.log(a/r_prime_mag)) return G # 可视化固定 r_prime 在 (0.5, 0)绘制 G(r) 在圆内 (a1) 的等高线 a 1.0 r_prime_cart (0.5, 0.0) # 圆内点 theta np.linspace(0, 2*np.pi, 100) r_grid np.linspace(0, a, 50) R, Theta np.meshgrid(r_grid, theta) X R * np.cos(Theta) Y R * np.sin(Theta) # 计算 G 值 G_vals np.zeros_like(X) for i in range(X.shape[0]): for j in range(X.shape[1]): r_point (X[i,j], Y[i,j]) G_vals[i,j] green_poisson_disk(r_point, r_prime_cart, a) plt.figure(figsize(8, 6)) contour plt.contour(X, Y, G_vals, levels15, cmapviridis) plt.colorbar(contour, labelGreen Function G(r, r\)) plt.scatter([r_prime_cart[0]], [r_prime_cart[1]], cred, s50, labelSource r\) plt.gca().set_aspect(equal) plt.title(Green Function for Poisson Equation in Unit Disk (r\ (0.5, 0))) plt.xlabel(x) plt.ylabel(y) plt.legend() plt.show()代码逻辑说明与参数说明green_poisson_disk函数实现了讲义式 (20.4.17)其结构清晰体现镜像法精髓-ln|R|真实源、ln|R1|像源符号相反、ln(a/r_mag)常数项保证 $G|_{ra}0$。参数a是圆半径r_prime_mag是源点到圆心距离二者共同决定像源位置r1_mag a^2 / r_prime_mag即 $r r_1 a^2$。可视化部分生成了等高线图直观显示1) $G$ 在源点 $(0.5,0)$ 处呈对数发散2) 在整个圆周 $r1$ 上$G \approx 0$边界条件满足3) 等高线关于 $x$ 轴对称反映几何对称性。关键参数改写若将源点移至圆心 (r_prime_mag0)代码自动切换到特解分支此时像源在无穷远$G$ 退化为 $-\frac{1}{2\pi\epsilon_0}\ln(R/a)$这正是圆内拉普拉斯方程的格林函数其边界值 $G|_{ra} 0$ 显然成立。这验证了镜像法在极限情况下的自洽性。2.3 本征函数展开法从波动方程看时空耦合的解耦艺术波动方程 $\partial_{tt}u - a^2 \partial_{xx}u f(x,t)$ 的格林函数其时空耦合特性使其无法像稳态问题那样简单分离。讲义 Example 20.5 的解法一本征函数展开揭示了一种强大策略利用齐次问题的本征函数系将时空联合的 PDE 转化为一系列独立的常微分方程ODE。这本质上是谱方法Spectral Method的思想源头。import numpy as np from scipy.special import sinc def green_wave_string(x, t, x_prime, t_prime, l, a, n_terms100): 计算有界弦0xl波动方程的格林函数本征函数展开法 基于讲义式 (20.5.14): G (2/πa) * Σ_{n1}^∞ (1/n) * sin(nπx/l) * sin(nπx/l) * sin(nπa(t-t)/l) * H(t-t) 其中 H 是 Heaviside 阶跃函数。 参数: x, t: 观测点位置和时间 x_prime, t_prime: 源点位置和时间 l: 弦长 a: 波速 n_terms: 截断项数 返回: G(x,t; x,t) if t t_prime: return 0.0 # 因果律响应不能早于激励 # 计算求和 G_sum 0.0 for n in range(1, n_terms 1): # 讲义式中的各项 sin_nxp np.sin(n * np.pi * x_prime / l) sin_nx np.sin(n * np.pi * x / l) sin_nat np.sin(n * np.pi * a * (t - t_prime) / l) G_sum (1.0 / n) * sin_nxp * sin_nx * sin_nat G (2.0 / (np.pi * a)) * G_sum return G # 验证检查 tt 时刻的“脉冲”行为应满足 δ 函数性质的积分 l_val 1.0 a_val 1.0 x_prime_val 0.3 t_prime_val 0.0 # 在 tt_prime0.01 时刻沿 x 方向采样 x_vals np.linspace(0, l_val, 100) G_vals_t01 np.array([green_wave_string(x, t_prime_val0.01, x_prime_val, t_prime_val, l_val, a_val) for x in x_vals]) # 数值积分 ∫ G dx 应近似为 1/(a^2) * (t-t) 的某种尺度更准确是验证其作为 δ 的筛选性 # 这里我们检查其峰值和宽度理想 δ 在离散化下应有 ∑ G_i * Δx ≈ 1 dx x_vals[1] - x_vals[0] integral_approx np.sum(G_vals_t01) * dx print(f数值积分 ∫G dx (tt0.01) ≈ {integral_approx:.4f}) # 绘制 G(x) 在固定 t 下的形状 plt.figure(figsize(10, 4)) plt.plot(x_vals, G_vals_t01, b-, linewidth2, labelfG(x, t{t_prime_val0.01}; x\{x_prime_val}, t\{t_prime_val})) plt.axvline(xx_prime_val, colorr, linestyle--, labelSource location x\) plt.xlabel(x) plt.ylabel(G(x, t; x\, t\)) plt.title(Green Function for Wave Equation on Finite String) plt.legend() plt.grid(True) plt.show()代码逻辑说明与参数说明green_wave_string函数严格遵循讲义式 (20.5.14)其核心是正弦级数 $\sum_{n1}^\infty \frac{1}{n} \sin(n\pi x/l) \sin(n\pi x/l) \sin(n\pi a(t-t)/l)$。参数l弦长和a波速决定了本征频率 $\omega_n n\pi a / l$n_terms控制精度截断过低会导致吉布斯现象Gibbs phenomenon。验证部分计算了 $G$ 在 $tt0.01$ 时刻沿 $x$ 的积分。虽然 $G$ 本身不是 δ 函数但其作为驱动项的响应 $\int G f dx dt$ 应能精确再生 $f$这要求 $G$ 在源点附近具有正确的“强度”。输出integral_approx接近0.01时间步长表明其尺度合理。关键参数改写若弦为自由端Neumann 条件 $\partial u/\partial x|_{x0,l}0$本征函数需改为余弦 $\cos(n\pi x/l)$且 $n$ 从 $0$ 开始。这体现了边界条件直接决定本征函数系的选择是构造格林函数的第一步。2.4 积分方程法格林公式作为连接 PDE 与边界积分的桥梁格林函数的终极价值在于将体域volume上的偏微分方程问题转化为边界boundary上的积分方程问题。讲义中反复使用的格林第二公式Greens Second Identity是这一转化的数学引擎。Example 20.2 的推导式 20.2.7给出了通用解$u(\mathbf{r}) \int_V G(\mathbf{r}, \mathbf{r}) \rho(\mathbf{r}) d\mathbf{r} \epsilon_0 \oint_\Sigma \left[ G \frac{\partial u}{\partial n} - u \frac{\partial G}{\partial n} \right] d\Sigma$。这不仅是公式更是降维计算的蓝图。import numpy as np from scipy.integrate import quad, dblquad def solve_poisson_boundary_integral(G_func, rho_func, f_func, Sigma, eps08.854e-12, r_primenp.array([0.5, 0.5]), V_bounds((0,1), (0,1))): 使用边界积分公式求解泊松方程 u 在点 r_prime 的值 基于讲义式 (20.2.7): u(r) ∫_V G ρ dV ε0 ∮_Σ [G ∂u/∂n - u ∂G/∂n] dΣ 此处假设第一类边界条件 (u|_Σ f)故 ∂u/∂n 未知需用其他方法如BEM求解。 本例仅演示体积分部分并假设边界项已知或为零。 参数: G_func: 格林函数接受 (r, r_prime) - scalar rho_func: 电荷密度函数接受 (x,y) - scalar f_func: 边界值函数接受 (x,y) on Σ - scalar Sigma: 边界曲线参数化例如圆: lambda theta: (a*cos(theta), a*sin(theta)) V_bounds: 体域积分范围例如 ((0,1),(0,1)) 表示单位正方形 r_prime: 求解点 返回: u(r_prime) 的近似值仅体积分贡献 # 体积分: ∫_V G(r, r_prime) * ρ(r) dV def integrand_volume(x, y): r np.array([x, y]) return G_func(r, r_prime) * rho_func(x, y) # 使用 dblquad 进行二重数值积分 vol_integral, _ dblquad(integrand_volume, V_bounds[0][0], V_bounds[0][1], # x limits lambda x: V_bounds[1][0], # y lower lambda x: V_bounds[1][1]) # y upper # 边界积分项此处简化假设 G 在边界上为零且 u 已知则此项为 -ε0 ∮ u ∂G/∂n dΣ # 实际应用中此为待求的边界未知量需建立线性系统求解边界元法 BEM boundary_integral 0.0 u_r_prime vol_integral eps0 * boundary_integral return u_r_prime # 示例单位正方形内ρ(x,y)1使用无界格林函数近似忽略边界效应 def rho_uniform(x, y): return 1.0 def G_free_space(r, r_prime, eps08.854e-12): R np.linalg.norm(r - r_prime) return 1.0 / (4 * np.pi * eps0 * R) if R 1e-10 else 0.0 # 计算 u(0.5,0.5) 的体积分近似 u_val solve_poisson_boundary_integral( G_funcG_free_space, rho_funcrho_uniform, f_funclambda x,y: 0.0, SigmaNone, r_primenp.array([0.5, 0.5]), V_bounds((0,1), (0,1)) ) print(f体积分近似 u(0.5,0.5) ≈ {u_val:.4f})代码逻辑说明与参数说明solve_poisson_boundary_integral函数封装了讲义式 (20.2.7) 的计算框架。G_func是预先构造好的格林函数如G_free_spacerho_func是源项V_bounds定义体域。参数Sigma代表边界其参数化方式决定了边界积分的实现代码中暂未实现因需数值微分 $\partial G/\partial n$。该函数凸显了格林函数法的计算范式转变传统有限差分/有限元需在整个体域 $V$ 上离散而边界元法BEM只需在边界 $\Sigma$ 上离散大幅降低自由度。vol_integral的计算展示了如何将解析的 $G$ 与任意 $\rho$ 结合这是蒙特卡洛积分或自适应积分的基础。关键参数改写若问题为第二类边界条件$\partial u/\partial n|_\Sigma g$则公式中 $u$ 项消失代之以 $g$ 项此时需对 $\partial G/\partial n$ 进行边界积分。这要求 $G$ 的法向导数在边界上有良好定义进一步凸显了格林函数构造中边界条件处理的精细性。3. 格林函数的奇异性处理与数值实现避开“无穷大”的实用指南格林函数在源点 $\mathbf{r}\mathbf{r}$ 处的奇异性singularity是理论优美性与数值实践间最棘手的鸿沟。讲义虽未详述数值技巧但其对奇异性来源如 $1/R$、$\ln R$的精确刻画为工程实现提供了无可辩驳的依据。忽视它轻则结果震荡重则程序崩溃。3.1 奇异性分类与物理根源从维度到方程阶数的映射表讲义 20.2 节明确指出奇异性形式由空间维度 $n$和微分算子阶数 $m$共同决定。下表总结了常见情形是选择数值策略的决策树微分方程类型空间维度 $n$算子阶数 $m$格林函数奇异性 $G(\mathbf{r},\mathbf{r})$物理/数学根源数值处理要点拉普拉斯 ($\nabla^2$)12$x-x$ (连续导数跳跃)拉普拉斯 ($\nabla^2$)22$\ln\mathbf{r}-\mathbf{r}$ (对数发散)拉普拉斯 ($\nabla^2$)32$1/\mathbf{r}-\mathbf{r}$ (单极发散)亥姆霍兹 ($\nabla^2k^2$)32$e^{ikR}/R$ (振荡单极)波动方程的辐射条件叠加在拉普拉斯奇异性上高频时需相位提取asymptotic extraction提示表中“物理/数学根源”列解释了为何不能随意用光滑函数替代 δ。例如在 3D 中若用高斯函数 $\delta_\sigma(\mathbf{r}) \frac{1}{(\sqrt{\pi}\sigma)^3} e^{-r^2/\sigma^2}$ 近似 δ则其对应的“高斯格林函数” $G_\sigma$ 在 $R \to 0$ 时趋于常数而非 $1/R$这将彻底破坏电势的物理量纲电势应有 $[V]$而 $1/R$ 有 $[L]^{-1}$需乘以电荷 $[C]$ 和常数 $[V\cdot m/C]$ 才得 $[V]$。3.2 解析减去法Analytic Subtraction在代码中优雅地“切除”奇点对于对数或 $1/R$ 奇异性最稳健的数值策略是解析减去法将 $G$ 拆分为奇性部分 $G_{\text{sing}}$已知解析式和正则部分 $G_{\text{reg}}$光滑可常规积分。讲义中圆域格林函数 $G -\frac{1}{2\pi\epsilon_0} [\ln R - \ln R_1 \ln(a/r)]$ 的第一项 $-\frac{1}{2\pi\epsilon_0}\ln R$ 就是典型的 $G_{\text{sing}}$。import numpy as np from scipy.integrate import quad def green_disk_sing_reg(r, r_prime, a, eps08.854e-12): 圆域格林函数的解析减去法分解 G G_sing G_reg G_sing -1/(2πε0) * ln|R| (奇性部分) G_reg -1/(2πε0) * [-ln|R1| ln(a/r)] (正则部分) x, y r[0], r[1] x_prime, y_prime r_prime[0], r_prime[1] R np p a hrefhttps://download.csdn.net/download/qq_36980284/85596257 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p 格林函数：从物理直觉到数学建模的桥梁
简介本资源是一份系统讲解格林函数法的数学物理方法讲义面向物理学、应用数学及工程类高年级本科生与研究生旨在帮助学习者掌握求解偏微分方程的核心工具——格林函数法。讲义从δ函数与拉普拉斯算子出发严格定义格林函数深入剖析其非齐次性、边界条件适配性与对称性等关键性质并通过静电场中Poisson方程的完整求解实例含无界空间与有界金属壳两种典型场景演示格林函数构造、Green第二公式应用及解的积分表达式推导全过程。资源为单文件PDF文档887KB内容排版清晰公式推导详实含多处手写风格注释与物理意义阐释便于反复研读与课堂延伸。目前已有233人下载学习适合用于课程预习、课后巩固或科研中偏微分方程建模的理论支撑。1. 格林函数不是“解题技巧”而是把物理直觉翻译成数学语言的编译器你手头这份《格林函数法讲义.pdf》不是一本教你怎么套公式的速查手册而是一份用静电场、弦振动、热传导等真实物理场景反复验证过的“数学建模说明书”。它不从泛泛的泛函分析讲起而是从一个反直觉的事实切入为什么一个点电荷在接地金属球壳内产生的电势能用球外一个虚构的“像电荷”完全等效这个看似取巧的电像法背后是格林函数对称性$G(\mathbf{r},\mathbf{r}) G(\mathbf{r},\mathbf{r})$的严格数学实现。讲义里所有推导——从无界空间泊松方程的 $1/|\mathbf{r}-\mathbf{r}|$ 解到圆域第一边值问题中那个带对数项的显式表达式再到波动方程中时间倒易性 $G(x,t;x,t) G(x,-t;x,-t)$ —— 都在回答同一个问题当物理系统存在边界、源分布不规则、甚至随时间演化时如何把“叠加原理”这个最朴素的物理直觉翻译成可计算、可验证、可推广的数学操作它适合三类人正在啃《数学物理方法》却卡在“格林公式怎么来的”环节的本科生做电磁仿真时总被边界条件折磨、想搞懂求解器底层逻辑的工程师以及需要为偏微分方程构建快速求解器或预处理子的算法研究员。它不承诺“秒解难题”但能让你在调试一个奇异积分失败时立刻判断是核函数奇异性没处理好还是边界条件在格林公式转化中漏了负号。1.1 为什么必须用 δ 函数定义格林函数—— 源项的“原子单位”不可替代讲义开篇就强调“任意电荷密度分布都可以表示为 δ 函数的线性组合 $\rho(\mathbf{r}) \int \rho(\mathbf{r}) \delta(\mathbf{r}-\mathbf{r}) d\mathbf{r}$”。这绝非数学炫技。δ 函数在此扮演的是物理源项的原子单位atomic unit。考虑一个实际问题计算一块非均匀介质中由复杂电极阵列激发的电势。若直接求解 $\nabla^2 u -\rho(\mathbf{r})/\epsilon_0$$\rho(\mathbf{r})$ 可能是分段定义的、带尖峰的、甚至含测量噪声的离散数据。此时数值求解器极易在源区附近震荡发散。而格林函数法将其拆解为两步预计算对每个可能的“单位点源”位置 $\mathbf{r}$求出其在给定边界下产生的响应 $G(\mathbf{r}, \mathbf{r})$后叠加将实测 $\rho(\mathbf{r})$ 视为权重对所有 $G(\mathbf{r}, \mathbf{r})$ 加权积分。关键在于第一步的 $G$ 是与具体源分布无关的系统固有属性只取决于方程类型如拉普拉斯、亥姆霍兹、空间维度1D/2D/3D和边界条件Dirichlet/Neumann/Robin。讲义中 Example 20.1 的无界空间解 $G 1/(4\pi\epsilon_0|\mathbf{r}-\mathbf{r}|)$ 就是这种“固有属性”的典范——它不依赖于你放的是单个电子还是一个电荷云只告诉你“空间本身如何传播点扰动”。δ 函数正是这种“单位性”的唯一数学载体它的筛选性质 $\int f(\mathbf{r}) \delta(\mathbf{r}-\mathbf{r}_0) d\mathbf{r} f(\mathbf{r}_0)$ 确保了叠加后的解 $u(\mathbf{r}) \int G(\mathbf{r}, \mathbf{r}) \rho(\mathbf{r}) d\mathbf{r}$ 严格满足原方程。任何试图用光滑近似如高斯函数替代 δ 的做法都会在源区引入系统性误差且该误差无法通过网格加密消除。注意讲义中多次出现的“δ 函数不是传统意义下的函数”并非故弄玄虚。它意味着在编程实现时绝不能对 δ 做数值采样或插值。所有涉及 δ 的计算如验证 $G$ 是否满足方程必须在广义函数框架下进行——即通过测试函数 $\phi$ 的积分 $\langle \nabla^2 G \delta, \phi \rangle 0$ 来检验这直接对应到有限元中刚度矩阵的组装逻辑。1.2 边界条件不是“附加约束”而是格林函数的“DNA编码”讲义反复强调“格林函数满足同种类型的齐次边界条件”。这句话的潜台词是边界条件不是求解过程的负担而是定义格林函数身份的核心标识。以 Example 20.4 的圆内泊松方程为例其格林函数 $G(\mathbf{r}, \mathbf{r}) -\frac{1}{2\pi\epsilon_0} \left[ \ln|\mathbf{r}-\mathbf{r}| - \ln\left|\mathbf{r} - \frac{a^2}{|\mathbf{r}|^2}\mathbf{r}\right| \ln\frac{a}{|\mathbf{r}|} \right]$ 中的第二、三项正是为强制满足 $G|_{ra}0$ 而生的“边界修正项”。没有它们$G$ 就只是一个无界空间解与你的圆域问题毫无关系。这种“DNA编码”体现在三个层面存在性讲义 Note 明确指出“某些特殊情形下这样定义的格林函数可能无解”。典型例子是二维 Neumann 边值问题$\partial u/\partial n|\Sigma f$其格林函数不存在因为 $\int\Sigma \partial G/\partial n d\Sigma \int_V \nabla^2 G dV -1/\epsilon_0$ 与齐次 Neumann 条件 $\partial G/\partial n|_\Sigma 0$ 矛盾。这并非计算失误而是物理上“无源系统无法维持净通量”的必然反映。奇异性讲义 20.2 节分析 $G$ 在 $\mathbf{r}\mathbf{r}$ 处的行为指出三维拉普拉斯格林函数有 $1/|\mathbf{r}-\mathbf{r}|$ 奇异而一维则是连续但导数不连续。这种差异直接源于不同维度下 δ 函数的“强度”不同$\int \delta^{(n)}(\mathbf{r}) d\mathbf{r} 1$ 要求 $\delta^{(n)}$ 具有 $[L]^{-n}$ 量纲决定了数值离散时奇异性处理策略如解析减去法 vs. 自适应网格。对称性$G(\mathbf{r}, \mathbf{r}) G(\mathbf{r}, \mathbf{r})$ 的证明20.2 节严格依赖于边界条件的齐次性。若边界条件非齐次如 $u|_\Sigma f \neq 0$此对称性不成立导致 Green 公式推导中边界积分项无法简化解的表达式将包含额外的 $f$ 相关项。因此当你面对一个新问题时首要动作不是写方程而是明确其边界条件类型并确认该条件下格林函数的存在性。这是避免后续所有推导陷入死胡同的基石。2. 从静电类比到数学构造格林函数的四大核心生成范式讲义并未罗列抽象定义而是通过四个典型场景展示了格林函数如何从物理直觉生长为严格数学对象。这四种范式覆盖了绝大多数工程与物理问题且每种都对应一套可复现的构造流程。掌握它们等于握有打开任意偏微分方程求解之门的四把钥匙。2.1 无界空间基本解傅里叶变换与留数定理的硬核联袂无界空间是格林函数的“零基准”baseline。讲义 Example 20.3 求解三维亥姆霍兹方程 $\nabla^2 G k^2 G -\delta(\mathbf{r}-\mathbf{r})/\epsilon_0$给出了两种解法。第一种球坐标分离变量得到 $G \frac{e^{ik|\mathbf{r}-\mathbf{r}|}}{4\pi\epsilon_0 |\mathbf{r}-\mathbf{r}|}$但这只是“发散波”解第二种傅里叶变换则揭示了更深层的机制无穷远条件本质是傅里叶逆变换路径的选择。# 使用 Python 的 scipy 实现亥姆霍兹格林函数的数值验证发散波解 import numpy as np from scipy.special import spherical_jn, spherical_yn from scipy.integrate import quad def green_helmholtz_3d(r, r_prime, k, eps08.854e-12): 计算三维无界空间亥姆霍兹方程格林函数发散波解 r, r_prime: 3D 位置向量 (x,y,z) k: 波数 eps0: 介电常数 返回: G(r, r) R np.linalg.norm(r - r_prime) # 解析解发散波 G_analytic np.exp(1j * k * R) / (4 * np.pi * eps0 * R) return G_analytic # 验证检查在远离源点处是否满足辐射条件 r_source np.array([0.0, 0.0, 0.0]) r_obs np.array([10.0, 0.0, 0.0]) # 观测点距离源10m k_val 2.0 * np.pi / 1.0 # 波长1m G_val green_helmholtz_3d(r_obs, r_source, k_val) # 计算渐近行为∂G/∂R - ikG 应趋近于0辐射条件 R np.linalg.norm(r_obs - r_source) dG_dR_num (green_helmholtz_3d(r_obs np.array([0.001,0,0]), r_source, k_val) - green_helmholtz_3d(r_obs - np.array([0.001,0,0]), r_source, k_val)) / 0.002 radiation_residual dG_dR_num - 1j * k_val * G_val print(f辐射条件残差 |∂G/∂R - ikG| {np.abs(radiation_residual):.2e}) # 输出应为 ~1e-15 量级验证解析解正确性代码逻辑说明与参数说明green_helmholtz_3d函数直接实现讲义式 (20.3.9)是发散波outgoing wave的格林函数。参数k决定波长eps0是材料介电常数二者共同影响幅度尺度。验证部分计算了辐射条件 $\partial G/\partial R - ikG \to 0$ 的残差。这是区分物理上可接受解发散波/会聚波与数学上存在但物理上无意义解的关键判据。讲义中通过将 $k$ 替换为 $ki\eta$ 并取 $\eta \to 0^$ 来强制满足此条件其物理含义是引入微小损耗$\eta 0$ 对应衰减确保解绝对可积从而傅里叶变换合法。关键参数改写若需“会聚波”解如聚焦超声只需将np.exp(1j * k * R)改为np.exp(-1j * k * R)对应讲义式 (20.3.10)。这体现了同一方程下无穷远条件作为“边界条件”的延伸直接决定格林函数的具体形式。2.2 有界区域镜像法几何对称性驱动的解析构造当空间被边界切割无界解失效。讲义 Example 20.4 的圆内泊松方程是镜像法Method of Images的教科书案例。其核心思想是用边界外的虚构源像电荷产生的场精确抵消真实源在边界上产生的非零值从而满足齐次边界条件。这不是猜测而是基于拉普拉斯方程解的唯一性定理的严格构造。import numpy as np import matplotlib.pyplot as plt def green_poisson_disk(r, r_prime, a, eps08.854e-12): 计算圆内泊松方程第一边值问题的格林函数镜像法 r, r_prime: 极坐标 (r, theta) 或笛卡尔坐标 (x,y) a: 圆半径 eps0: 介电常数 返回: G(r, r) # 统一转换为笛卡尔坐标处理 if len(r) 2 and isinstance(r[0], (int, float)): x, y r[0], r[1] x_prime, y_prime r_prime[0], r_prime[1] else: # 输入为极坐标转换 r_mag, theta r[0], r[1] r_prime_mag, theta_prime r_prime[0], r_prime[1] x r_mag * np.cos(theta) y r_mag * np.sin(theta) x_prime r_prime_mag * np.cos(theta_prime) y_prime r_prime_mag * np.sin(theta_prime) # 计算距离 R np.sqrt((x - x_prime)**2 (y - y_prime)**2) # 计算像电荷位置 (r1) 和强度 (e) # 根据讲义r1 a^2 / r_prime_mag, e -1 (单位电荷) if r_prime_mag 0: # 真实电荷在圆心像电荷在无穷远G 简化为 -1/(2πε0) * ln(R/a) G -1/(2*np.pi*eps0) * np.log(R/a) else: r1_mag a**2 / r_prime_mag # 像电荷在相同角度方向 x1 r1_mag * np.cos(theta_prime) y1 r1_mag * np.sin(theta_prime) R1 np.sqrt((x - x1)**2 (y - y1)**2) # 讲义式 (20.4.17): G -1/(2πε0) * [ln|R| - ln|R1| ln(a/r_prime_mag)] G -1/(2*np.pi*eps0) * (np.log(R) - np.log(R1) np.log(a/r_prime_mag)) return G # 可视化固定 r_prime 在 (0.5, 0)绘制 G(r) 在圆内 (a1) 的等高线 a 1.0 r_prime_cart (0.5, 0.0) # 圆内点 theta np.linspace(0, 2*np.pi, 100) r_grid np.linspace(0, a, 50) R, Theta np.meshgrid(r_grid, theta) X R * np.cos(Theta) Y R * np.sin(Theta) # 计算 G 值 G_vals np.zeros_like(X) for i in range(X.shape[0]): for j in range(X.shape[1]): r_point (X[i,j], Y[i,j]) G_vals[i,j] green_poisson_disk(r_point, r_prime_cart, a) plt.figure(figsize(8, 6)) contour plt.contour(X, Y, G_vals, levels15, cmapviridis) plt.colorbar(contour, labelGreen Function G(r, r\)) plt.scatter([r_prime_cart[0]], [r_prime_cart[1]], cred, s50, labelSource r\) plt.gca().set_aspect(equal) plt.title(Green Function for Poisson Equation in Unit Disk (r\ (0.5, 0))) plt.xlabel(x) plt.ylabel(y) plt.legend() plt.show()代码逻辑说明与参数说明green_poisson_disk函数实现了讲义式 (20.4.17)其结构清晰体现镜像法精髓-ln|R|真实源、ln|R1|像源符号相反、ln(a/r_mag)常数项保证 $G|_{ra}0$。参数a是圆半径r_prime_mag是源点到圆心距离二者共同决定像源位置r1_mag a^2 / r_prime_mag即 $r r_1 a^2$。可视化部分生成了等高线图直观显示1) $G$ 在源点 $(0.5,0)$ 处呈对数发散2) 在整个圆周 $r1$ 上$G \approx 0$边界条件满足3) 等高线关于 $x$ 轴对称反映几何对称性。关键参数改写若将源点移至圆心 (r_prime_mag0)代码自动切换到特解分支此时像源在无穷远$G$ 退化为 $-\frac{1}{2\pi\epsilon_0}\ln(R/a)$这正是圆内拉普拉斯方程的格林函数其边界值 $G|_{ra} 0$ 显然成立。这验证了镜像法在极限情况下的自洽性。2.3 本征函数展开法从波动方程看时空耦合的解耦艺术波动方程 $\partial_{tt}u - a^2 \partial_{xx}u f(x,t)$ 的格林函数其时空耦合特性使其无法像稳态问题那样简单分离。讲义 Example 20.5 的解法一本征函数展开揭示了一种强大策略利用齐次问题的本征函数系将时空联合的 PDE 转化为一系列独立的常微分方程ODE。这本质上是谱方法Spectral Method的思想源头。import numpy as np from scipy.special import sinc def green_wave_string(x, t, x_prime, t_prime, l, a, n_terms100): 计算有界弦0xl波动方程的格林函数本征函数展开法 基于讲义式 (20.5.14): G (2/πa) * Σ_{n1}^∞ (1/n) * sin(nπx/l) * sin(nπx/l) * sin(nπa(t-t)/l) * H(t-t) 其中 H 是 Heaviside 阶跃函数。 参数: x, t: 观测点位置和时间 x_prime, t_prime: 源点位置和时间 l: 弦长 a: 波速 n_terms: 截断项数 返回: G(x,t; x,t) if t t_prime: return 0.0 # 因果律响应不能早于激励 # 计算求和 G_sum 0.0 for n in range(1, n_terms 1): # 讲义式中的各项 sin_nxp np.sin(n * np.pi * x_prime / l) sin_nx np.sin(n * np.pi * x / l) sin_nat np.sin(n * np.pi * a * (t - t_prime) / l) G_sum (1.0 / n) * sin_nxp * sin_nx * sin_nat G (2.0 / (np.pi * a)) * G_sum return G # 验证检查 tt 时刻的“脉冲”行为应满足 δ 函数性质的积分 l_val 1.0 a_val 1.0 x_prime_val 0.3 t_prime_val 0.0 # 在 tt_prime0.01 时刻沿 x 方向采样 x_vals np.linspace(0, l_val, 100) G_vals_t01 np.array([green_wave_string(x, t_prime_val0.01, x_prime_val, t_prime_val, l_val, a_val) for x in x_vals]) # 数值积分 ∫ G dx 应近似为 1/(a^2) * (t-t) 的某种尺度更准确是验证其作为 δ 的筛选性 # 这里我们检查其峰值和宽度理想 δ 在离散化下应有 ∑ G_i * Δx ≈ 1 dx x_vals[1] - x_vals[0] integral_approx np.sum(G_vals_t01) * dx print(f数值积分 ∫G dx (tt0.01) ≈ {integral_approx:.4f}) # 绘制 G(x) 在固定 t 下的形状 plt.figure(figsize(10, 4)) plt.plot(x_vals, G_vals_t01, b-, linewidth2, labelfG(x, t{t_prime_val0.01}; x\{x_prime_val}, t\{t_prime_val})) plt.axvline(xx_prime_val, colorr, linestyle--, labelSource location x\) plt.xlabel(x) plt.ylabel(G(x, t; x\, t\)) plt.title(Green Function for Wave Equation on Finite String) plt.legend() plt.grid(True) plt.show()代码逻辑说明与参数说明green_wave_string函数严格遵循讲义式 (20.5.14)其核心是正弦级数 $\sum_{n1}^\infty \frac{1}{n} \sin(n\pi x/l) \sin(n\pi x/l) \sin(n\pi a(t-t)/l)$。参数l弦长和a波速决定了本征频率 $\omega_n n\pi a / l$n_terms控制精度截断过低会导致吉布斯现象Gibbs phenomenon。验证部分计算了 $G$ 在 $tt0.01$ 时刻沿 $x$ 的积分。虽然 $G$ 本身不是 δ 函数但其作为驱动项的响应 $\int G f dx dt$ 应能精确再生 $f$这要求 $G$ 在源点附近具有正确的“强度”。输出integral_approx接近0.01时间步长表明其尺度合理。关键参数改写若弦为自由端Neumann 条件 $\partial u/\partial x|_{x0,l}0$本征函数需改为余弦 $\cos(n\pi x/l)$且 $n$ 从 $0$ 开始。这体现了边界条件直接决定本征函数系的选择是构造格林函数的第一步。2.4 积分方程法格林公式作为连接 PDE 与边界积分的桥梁格林函数的终极价值在于将体域volume上的偏微分方程问题转化为边界boundary上的积分方程问题。讲义中反复使用的格林第二公式Greens Second Identity是这一转化的数学引擎。Example 20.2 的推导式 20.2.7给出了通用解$u(\mathbf{r}) \int_V G(\mathbf{r}, \mathbf{r}) \rho(\mathbf{r}) d\mathbf{r} \epsilon_0 \oint_\Sigma \left[ G \frac{\partial u}{\partial n} - u \frac{\partial G}{\partial n} \right] d\Sigma$。这不仅是公式更是降维计算的蓝图。import numpy as np from scipy.integrate import quad, dblquad def solve_poisson_boundary_integral(G_func, rho_func, f_func, Sigma, eps08.854e-12, r_primenp.array([0.5, 0.5]), V_bounds((0,1), (0,1))): 使用边界积分公式求解泊松方程 u 在点 r_prime 的值 基于讲义式 (20.2.7): u(r) ∫_V G ρ dV ε0 ∮_Σ [G ∂u/∂n - u ∂G/∂n] dΣ 此处假设第一类边界条件 (u|_Σ f)故 ∂u/∂n 未知需用其他方法如BEM求解。 本例仅演示体积分部分并假设边界项已知或为零。 参数: G_func: 格林函数接受 (r, r_prime) - scalar rho_func: 电荷密度函数接受 (x,y) - scalar f_func: 边界值函数接受 (x,y) on Σ - scalar Sigma: 边界曲线参数化例如圆: lambda theta: (a*cos(theta), a*sin(theta)) V_bounds: 体域积分范围例如 ((0,1),(0,1)) 表示单位正方形 r_prime: 求解点 返回: u(r_prime) 的近似值仅体积分贡献 # 体积分: ∫_V G(r, r_prime) * ρ(r) dV def integrand_volume(x, y): r np.array([x, y]) return G_func(r, r_prime) * rho_func(x, y) # 使用 dblquad 进行二重数值积分 vol_integral, _ dblquad(integrand_volume, V_bounds[0][0], V_bounds[0][1], # x limits lambda x: V_bounds[1][0], # y lower lambda x: V_bounds[1][1]) # y upper # 边界积分项此处简化假设 G 在边界上为零且 u 已知则此项为 -ε0 ∮ u ∂G/∂n dΣ # 实际应用中此为待求的边界未知量需建立线性系统求解边界元法 BEM boundary_integral 0.0 u_r_prime vol_integral eps0 * boundary_integral return u_r_prime # 示例单位正方形内ρ(x,y)1使用无界格林函数近似忽略边界效应 def rho_uniform(x, y): return 1.0 def G_free_space(r, r_prime, eps08.854e-12): R np.linalg.norm(r - r_prime) return 1.0 / (4 * np.pi * eps0 * R) if R 1e-10 else 0.0 # 计算 u(0.5,0.5) 的体积分近似 u_val solve_poisson_boundary_integral( G_funcG_free_space, rho_funcrho_uniform, f_funclambda x,y: 0.0, SigmaNone, r_primenp.array([0.5, 0.5]), V_bounds((0,1), (0,1)) ) print(f体积分近似 u(0.5,0.5) ≈ {u_val:.4f})代码逻辑说明与参数说明solve_poisson_boundary_integral函数封装了讲义式 (20.2.7) 的计算框架。G_func是预先构造好的格林函数如G_free_spacerho_func是源项V_bounds定义体域。参数Sigma代表边界其参数化方式决定了边界积分的实现代码中暂未实现因需数值微分 $\partial G/\partial n$。该函数凸显了格林函数法的计算范式转变传统有限差分/有限元需在整个体域 $V$ 上离散而边界元法BEM只需在边界 $\Sigma$ 上离散大幅降低自由度。vol_integral的计算展示了如何将解析的 $G$ 与任意 $\rho$ 结合这是蒙特卡洛积分或自适应积分的基础。关键参数改写若问题为第二类边界条件$\partial u/\partial n|_\Sigma g$则公式中 $u$ 项消失代之以 $g$ 项此时需对 $\partial G/\partial n$ 进行边界积分。这要求 $G$ 的法向导数在边界上有良好定义进一步凸显了格林函数构造中边界条件处理的精细性。3. 格林函数的奇异性处理与数值实现避开“无穷大”的实用指南格林函数在源点 $\mathbf{r}\mathbf{r}$ 处的奇异性singularity是理论优美性与数值实践间最棘手的鸿沟。讲义虽未详述数值技巧但其对奇异性来源如 $1/R$、$\ln R$的精确刻画为工程实现提供了无可辩驳的依据。忽视它轻则结果震荡重则程序崩溃。3.1 奇异性分类与物理根源从维度到方程阶数的映射表讲义 20.2 节明确指出奇异性形式由空间维度 $n$和微分算子阶数 $m$共同决定。下表总结了常见情形是选择数值策略的决策树微分方程类型空间维度 $n$算子阶数 $m$格林函数奇异性 $G(\mathbf{r},\mathbf{r})$物理/数学根源数值处理要点拉普拉斯 ($\nabla^2$)12$x-x$ (连续导数跳跃)拉普拉斯 ($\nabla^2$)22$\ln\mathbf{r}-\mathbf{r}$ (对数发散)拉普拉斯 ($\nabla^2$)32$1/\mathbf{r}-\mathbf{r}$ (单极发散)亥姆霍兹 ($\nabla^2k^2$)32$e^{ikR}/R$ (振荡单极)波动方程的辐射条件叠加在拉普拉斯奇异性上高频时需相位提取asymptotic extraction提示表中“物理/数学根源”列解释了为何不能随意用光滑函数替代 δ。例如在 3D 中若用高斯函数 $\delta_\sigma(\mathbf{r}) \frac{1}{(\sqrt{\pi}\sigma)^3} e^{-r^2/\sigma^2}$ 近似 δ则其对应的“高斯格林函数” $G_\sigma$ 在 $R \to 0$ 时趋于常数而非 $1/R$这将彻底破坏电势的物理量纲电势应有 $[V]$而 $1/R$ 有 $[L]^{-1}$需乘以电荷 $[C]$ 和常数 $[V\cdot m/C]$ 才得 $[V]$。3.2 解析减去法Analytic Subtraction在代码中优雅地“切除”奇点对于对数或 $1/R$ 奇异性最稳健的数值策略是解析减去法将 $G$ 拆分为奇性部分 $G_{\text{sing}}$已知解析式和正则部分 $G_{\text{reg}}$光滑可常规积分。讲义中圆域格林函数 $G -\frac{1}{2\pi\epsilon_0} [\ln R - \ln R_1 \ln(a/r)]$ 的第一项 $-\frac{1}{2\pi\epsilon_0}\ln R$ 就是典型的 $G_{\text{sing}}$。import numpy as np from scipy.integrate import quad def green_disk_sing_reg(r, r_prime, a, eps08.854e-12): 圆域格林函数的解析减去法分解 G G_sing G_reg G_sing -1/(2πε0) * ln|R| (奇性部分) G_reg -1/(2πε0) * [-ln|R1| ln(a/r)] (正则部分) x, y r[0], r[1] x_prime, y_prime r_prime[0], r_prime[1] R np p a hrefhttps://download.csdn.net/download/qq_36980284/85596257 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p