用C语言求解像管电场：从拉普拉斯方程到连续超松弛迭代

简介：北京理工大学光电学院《光电成像器件计算机辅助设计CAD》课程实验报告，核心内容是基于电场高斯定律、环路定律与电势差定理，配合连续超张弛迭代算法与拉格朗日插值算法，用C语言实现像管内电势计算的CAD程序设计。报告完整记录了从网格划分、有限差分方程建立到超张弛迭代收敛、边界封闭处理的数值求解流程，并给出了圆柱形像管内任意点电位计算、对称轴电位拉格朗日插值的实现方法。资源共1个PDF文档，压缩包743KB，适合光电工程、微光夜视技术方向的高年级本科生与研究生阅读，也可供从事电子光学系统仿真及数值计算工作的工程技术人员参考。已有171人学习下载。通过这份实验报告，读者可掌握将电磁场理论模型转化为可运行程序的方法，理解CAD软件在像管电场分布计算与等势线绘制中的具体应用，为像管设计、电子轨迹分析与像质评估打下工程实践基础。

1. 像管电场计算的一种落地解法：用C语言把拉普拉斯方程变成可迭代的CAD工具

拿到一份像管设计图，阴极、阳极、荧光屏，电极形状不规则，手写拉普拉斯方程去求解析解几乎是不可能完成的任务。解析法只对少数对称理想模型有效，实验测量法又受制于探头尺寸和场扰动，精度很难上去。这时候，用C语言写一个数值求解器：在网格上把偏微分方程离散成差分方程，通过连续超松弛迭代在毫秒级逼近真实电场分布，再用拉格朗日插值求任意位置的电势，就是一条非常务实的工程路径。这个项目正是一个完整的最小实现，代码量不算大，却把有限差分、连续超松弛迭代、等势线扫描绘制、叠加定理验证串成了一条完整的计算机辅助设计工具链。它不仅算出了像管内部的电势分布，还把迭代因子、残差、等势线坐标全部落盘到文本文件，非常适合想搞懂“电场数值计算到底怎么落地成软件”的开发者，尤其是对电磁场数值解法刚上手的研究生或工程师。

2. 理论基础与整体架构：从三大静电场定理到五段式模块设计

2.1 静电场三大定理为什么是数值算法的“免死金牌”

在着手写代码之前，必须先把物理边界条件想清楚。静电场有三大基本定理，这不仅是理论推导的基石，更是你调试程序时判断结果对错的标尺。

第一个是唯一性定理。如果给定系统中所有电极的形状、排列以及每个电极的电位，那么由这些电极产生的静电场将由拉普拉斯方程唯一确定。这句话的工程含义很直接：只要你的网格足够密、边界条件给得正确，不管你用迭代法还是直接法求解，得到的电位分布应该收敛到同一个值。如果你在程序中尝试了多种初始猜测值，最后迭代出来的结果不一样，那一定是你的边界封闭处理或者差分格式写错了，而不是算法发散的问题。

第二个是相似性定理。若各电极的电位都增大为K倍，空间各点的电位也增大为K倍，等势面形状不变。这个定理在验证程序时非常有用：你把输入文件里所有电极电压都乘2，算出来的每个网格点电压也应该乘2，等势线坐标则应该完全相同。如果做不到这一点，说明你的程序中存在某些“硬编码”的常数影响了结果的线性比例关系。

第三个是多电极系统的电位叠加定理。当各电极形状、相对位置确定后，总电场分布等于各电极单独作用（其余电极接地）时电场的线性叠加。这个定理是程序中“叠加定理验证模块”的核心依据。实现上就是跑n遍单电极的迭代，把结果加起来，再跟一次跑完所有电极的完整结果做逐点对比。这既是对算法的验证，也是对边界封闭处理正确性的交叉检查。

2.2 不等距五点差分格式：从拉普拉斯方程到线性方程组

在无空间电荷的真空中，静电场电位满足拉普拉斯方程。解析求解时，这是偏微分方程；但数值求解时，需要把它离散成差分方程。程序采用旋转对称场下的“十”字形不等距五点差分格式。假设网格点坐标为(z, r)，电位为φ0，相邻上下左右四个点的间距分别为h1、h2、h3、h4（取绝对值），则差分公式可以写成：

φ0 = (c1*φ1 + c2*φ2 + c3*φ3 + c4*φ4) / c0

其中系数c1到c4由相邻点的间距决定，c0是它们的加权和。对于轴上点（r=0），由于径向方向的间距存在0点，需要用洛必达法则处理c3、c4的计算公式。这部分是程序里最容易写错的地方，因为一旦网格划分不均匀，系数的计算就不能用简单的“等间距差分模板”套，必须根据实际网格坐标动态计算各点的间距和权重。

实现时我的习惯是单独写一个系数计算函数，输入是当前网格点的四个邻居坐标，输出是c0到c4。这样等间距网格和不等距网格可以共用同一套差分迭代逻辑，后续如果要换网格划分方式，就不需要动主迭代循环了。

2.3 程序模块划分：输入、初始化、迭代、绘图、验证五段式

这个程序的整体架构是典型的“串行流水线”风格，模块边界非常清晰，非常像早期C语言工程的标准组织方式。全程序一共分成五个功能模块：

• 输入模块：负责读入像管参数、初始化部分变量，并原样输出参数到日志文件。
• 初始化模块：负责划分网格步长，初始化电位矩阵，用线性插值法封闭像管边界。
• 超松弛迭代模块：这是整个CAD程序的心脏，负责对像管内部电势分布进行数值计算。
• 绘图与输出模块：进行行列扫描得到等位线坐标，用EasyX插件绘制电极和等位线图形，再进行拉格朗日插值绘制轴上电位变化曲线。
• 叠加定理验证模块：将每个电极单独置为有效电位，其余电极为0，调用迭代模块计算单位电位分布函数，最后叠加并与完整计算结果对比。

这种五段式架构的优势在于“可替换性”。比如你不想用EasyX绘图，只需要改绘图模块，其他模块完全不用动；你想把迭代算法换成共轭梯度法，也只需要替换超松弛迭代模块本身，输入输出格式和行为完全兼容。

3. 核心代码拆解：SOR迭代器与等势线扫描的实现细节

3.1 差分方程系数计算与边界封闭的C语言实现

这里给出系数计算的简化版本。需要注意的是，实际工程代码通常会预先计算好所有网格点上的系数，并存入数组中，避免每次迭代都重复计算，因为这是典型的“内存换时间”优化策略。

// 计算不等距五点差分格式的系数 // 输入：当前网格点坐标 (z, r) 以及上下左右四个邻居的坐标 // 输出：差分系数 cc[0] 到 cc[4]，其中 cc[0] 是中心点权重分母 void calc_coeff(double z, double r, double h1, double h2, double h3, double h4, double *cc) { double c1, c2, c3, c4; // 轴向差分系数：直接按等距公式计算 c1 = 2.0 / (h1 * (h1 + h2)); c2 = 2.0 / (h2 * (h1 + h2)); // 径向差分系数：需要处理 r=0 的奇异性 if (r < 1e-9) { // 轴上点，使用洛必达法则后的极限公式 c3 = 4.0 / (h3 * (h3 + h4)); c4 = 4.0 / (h4 * (h3 + h4)); } else { // 非轴上点，常规公式 c3 = (2.0 * r + h3) / (h3 * r * (h3 + h4)); c4 = (2.0 * r - h4) / (h4 * r * (h3 + h4)); } cc[0] = c1 + c2 + c3 + c4; cc[1] = c1; cc[2] = c2; cc[3] = c3; cc[4] = c4; }

代码里的关键点在于径向系数的来源。在柱坐标下拉普拉斯算子的径向部分带有1/r项，直接做差分会在r=0处产生奇点。程序的处理方式是对轴上点单独推导极限公式，这能保证对称轴上的电位计算精度不会因为坐标奇异性而明显劣化。实际调试时，如果发现轴上电位曲线出现毛刺或不光滑，先检查这一段的系数计算是否正确。

边界封闭的处理同样重要。像管电极并不是完全封闭的，电极之间存在空白区域，这些区域的边界电位是未知的。如果不做处理，迭代计算就缺少边界值，无法收敛。程序采用最简单的线性插值法：在两个相邻电极的边界之间，按照距离比例给空白区域内的边界点赋予电位值。这样做带来的误差在实际像管设计中是可以接受的，因为远离电极的边界区域对轴上电位的影响本身就很小。

3.2 连续超松弛迭代核心循环：从初始猜测到动态调整

连续超松弛迭代的核心思想是：每次遍历网格时，用当前点的四个邻居更新值来计算新的电位，再引入一个松弛因子ω来加速收敛。当1 < ω < 2时，收敛过程会明显加快，但ω过大会导致震荡甚至发散。

程序中的迭代策略分为三个阶段：第一阶段用ω=1进行纯高斯-赛德尔迭代，作为基准启动；第二阶段用固定的ω=1.375迭代12次；之后每12次迭代作为一个“轮次”，在轮次结束时重新计算最优松弛因子，并更新全局的ω值。

// 连续超松弛迭代核心循环 // phi: 电位矩阵, R: 行数, C: 列数 // eps: 迭代控制精度, max_round: 最大轮次 void sor_solve(double **phi, int R, int C, double eps, int max_round) { double omega = 1.0; double max_res = 0.0, avg_res = 0.0; for (int rnd = 1; rnd <= max_round; rnd++) { // 每轮迭代次数，前两轮特殊处理，后续每轮12次 int iters_per_round = (rnd <= 2) ? (rnd == 1 ? 1 : 12) : 12; for (int k = 0; k < iters_per_round; k++) { max_res = 0.0; avg_res = 0.0; int cnt = 0; // 遍历所有内部网格点 for (int i = 1; i < R - 1; i++) { for (int j = 1; j < C - 1; j++) { double old = phi[i][j]; double new_val = 0.0; // 这里需要根据网格点类型选取对应的差分系数 double c0, c1, c2, c3, c4; get_coeff_at(i, j, &c0, &c1, &c2, &c3, &c4); new_val = (c1 * phi[i-1][j] + c2 * phi[i+1][j] + c3 * phi[i][j-1] + c4 * phi[i][j+1]) / c0; // SOR 更新公式 phi[i][j] = (1.0 - omega) * old + omega * new_val; double res = fabs(phi[i][j] - old); max_res = fmax(max_res, res); avg_res += res; cnt++; } } avg_res /= cnt; // 输出本轮迭代的残差 printf("第 %d 轮，第 %d 次迭代，迭代因子: %f 平均残差：%f , 最大残差：%f\n", rnd, k + 1, omega, avg_res, max_res); // 满足精度要求则提前退出 if (max_res < eps) { printf("迭代收敛！\n"); return; } } // 每轮迭代结束后更新松弛因子 // 实际代码中omega的计算会用到谱半径估计，这里简化为递增策略 if (rnd == 1) omega = 1.375; else if (rnd == 2) omega = 1.75; else omega = fmin(1.95, omega + 0.05); } }

这段代码最需要注意的地方有两个。第一，new_val的更新必须使用“最新值”，也就是当前轮次已经计算过的节点值，这是高斯-赛德尔迭代的基本要求。如果你用上一轮次的旧值做更新，那就变成了雅可比迭代，收敛速度会大幅下降。第二，omega的更新策略需要根据实际残差变化做调整。在原始实验的迭代日志中可以看到，omega从1.0逐步爬到1.85左右，整个过程非常平滑。如果你在调试自己代码时发现omega调整后残差反而变大，说明当前问题域需要用更保守的更新策略。

3.3 行列扫描绘制等势线：把电位矩阵变成坐标点集

迭代完成之后，你手里有一整张网格的电位值。要画等势线，需要把“哪几个点电位相等”这个问题转化为算法可执行的操作。程序采用的方法非常经典：行扫描加列扫描，然后对相邻网格点做线性插值。

// 扫描并记录所有等势线点 // phi: 电位矩阵, V_level: 当前要搜索的电位值 // points: 输出点集, max_points: 点集容量 int scan_equipotential(double **phi, int R, int C, double V_level, Point *points, int max_points) { int count = 0; // 行扫描：沿着每一行检查水平方向穿越 for (int i = 1; i < R - 1; i++) { for (int j = 1; j < C - 1; j++) { if (phi[i][j] < V_level && phi[i][j+1] > V_level) { // 在线性插值求精确z坐标 double t = (V_level - phi[i][j]) / (phi[i][j+1] - phi[i][j]); double z = j + t * 1.0; // 假设网格间距为1 double r = i * 1.0; if (count < max_points) { points[count].z = z; points[count].r = r; count++; } } } } // 列扫描：沿着每一列检查垂直方向穿越 for (int j = 1; j < C - 1; j++) { for (int i = 1; i < R - 1; i++) { if (phi[i][j] < V_level && phi[i+1][j] > V_level) { // 在线性插值求精确r坐标 double t = (V_level - phi[i][j]) / (phi[i+1][j] - phi[i][j]); double r = i + t * 1.0; double z = j * 1.0; if (count < max_points) { points[count].z = z; points[count].r = r; count++; } } } } // 对点集按z坐标排序 qsort(points, count, sizeof(Point), compare_z); return count; }

行扫描和列扫描都会产生等势线点，但单独用其中一种会漏掉那些与网格线平行的等势线段。行列结合能最大程度保证等势线的连续性。扫描结束后得到的点集并不是有序的绘图路径，所以程序还需要对点集按照z坐标进行排序，然后按顺序连接。这个排序步骤在原始报告中描述得很清楚：先按照z坐标增大顺序排序，连接一次，再按照z坐标减小顺序排序，连接一次，主要目的是处理等势线在像管底部发生回折的情况。

4. 把程序跑起来：dat输入文件与EasyX绘图环境配置

4.1 输入文件参数详解：从电极个数到迭代精度

输入文件是后缀为.dat的文本文件，每一行对应一类参数。文件的组织方式几乎是“无格式”的纯序列文本，参数的顺序完全固定，任何一行读错了，后面的数据全部错位。

从输入格式可以看到，程序对网格划分的约束是M1 + M2 <= 50，所有电极厚度相同，所有电极内孔径半径相同。这些约束在原始代码里被写成了硬编码的数组大小，如果直接跑实验数据没问题，但想改成像管结构就需要动全局常量。

4.2 开发环境配置与编译运行步骤

程序使用C语言编写，绘图依赖EasyX插件，所以运行的推荐环境是Windows系统加Visual Studio。具体配置步骤为：第一步，下载并安装Visual Studio社区版；第二步，从官网下载EasyX库，运行安装程序，它会自动识别当前环境中的Visual Studio版本并完成集成；第三步，新建一个空项目，把源码文件加入项目，在项目属性中确认“字符集”设置为“使用多字节字符集”，因为早期代码习惯用char*处理文件名，如果是Unicode字符集可能会报类型不匹配的编译错误；第四步，把输入文件放到与main.c相同的目录下，直接编译运行。

运行时程序会在控制台窗口打印迭代过程，包括每一轮次、每一次迭代的迭代因子、平均残差和最大残差，就像附录一展示的那样。全部迭代完成后，会弹出EasyX绘图窗口显示等势线和轴上电位曲线，同时会在输入文件同目录下生成一个.txt文件，里面包含了完整的迭代结果和等势线坐标。由于绘图窗口是图形界面，程序退出后窗口会消失，如果看不到绘图结果，需要在main函数末尾加上system("pause");或者getchar();来保留窗口。

4.3 避坑指南：迭代不收敛、边界封闭失效与绘图黑屏

现象一：迭代到第3轮时残差非但不下降，反而突然增大。

原因：松弛因子ω更新过快，跳过了最优点，导致迭代过程进入震荡区间。原始日志中ω从1.375跳到1.755时，最大残差从0.103反弹到0.280，这是典型的过松弛震荡。解决：把ω的更新步长调小，每次只增加0.02到0.03，或者引入自适应机制，当检测到最大残差连续两次上升时，自动回退ω值。

现象二：绘制的等势线在电极之间的缺口处断开，或者明显扭曲。

原因：边界封闭处理不正确。如果两个相邻电极的电位差很大，线性插值得到的边界点电位可能不符合实际的物理衰减规律，导致靠近边界的网格点电位被错误牵引。解决：检查边界封闭的赋值顺序，确保左右边界、上下边界各自独立处理，不要让插值公式跨边界引用。另外可以尝试在电极之间增加少量过渡网格点，让插值更平缓。

现象三：EasyX绘图窗口一闪而过，根本来不及看结果。

原因：程序运行完main函数后，进程退出，图形窗口被系统自动回收。解决：在main函数的绘图代码后增加getchar();或system("pause");，让程序挂起等待键盘输入。如果是在Visual Studio里按F5运行，也可以改用Ctrl+F5（开始执行但不调试），这样控制台窗口结束后不会自动关闭。

现象四：修改网格数量M1或M2后，程序直接崩溃或抛出越界错误。

原因：原始代码的数组大小是按50x113硬编码的，M1+M2超过50时就会越界写入内存。解决：建议把这些数组定义改成动态分配，使用malloc和free，或者直接定义足够大的常量（比如200），并加上输入参数校验，在读取M1+M2时如果超过上限就报错退出，而不是等到运行到一半才崩溃。

5. 数据处理与误差分析：残差日志、叠加定理验证与实际工程意义

5.1 迭代残差日志解读：平均残差与最大残差的玄学

迭代收敛的过程不是一条光滑下降的曲线，而是带着“毛刺”的震荡衰减。查看附录一里的残差日志，你会发现第3轮第1次迭代的平均残差是0.056243，比第2轮最后一次迭代的0.021593还要大。这正是调整松弛因子后的“回弹”现象。松弛因子的调整会让算法暂时偏离原来的路径，但随后会以更快的速度收敛。

解读残差日志时，重点关注两点。第一是“最大残差”，它代表当前网格中最不收敛的那个点的误差大小。如果最大残差一直不降，说明某个区域可能被错误的边界条件污染了，比如边界封闭线上有某个点的电位值没被正确初始化，会持续向周围网格“发射”错误信号。第二是“平均残差”与“最大残差”的比值。正常情况下，这个比值应该在0.2到0.5之间。如果平均残差很小但最大残差很大，说明整个场分布已经基本稳定，但有几个孤立点存在异常，这种时候通常要回溯检查那几个点的坐标位置，看看是否落在了电极边缘或边界线上。

5.2 叠加定理验证模块：把n次SOR结果叠起来

叠加定理的验证逻辑非常朴素，但工程价值很高。具体实现是：先将所有电极电位清零，只保留第1个电极的电位，运行一次完整的连续超松弛迭代模块，把结果存入一个独立的叠加数组；然后清零所有电极，只保留第2个电极的电位，再运行一次迭代，把结果累加到刚才的数组里；重复这个过程，直到所有电极都单独处理完成。

// 叠加定理验证模块的核心伪代码 // base_phi: 正常计算得到的总电场电位 // super_phi: 叠加定理计算得到的电位 // single_phi: 单电极电场电位数组 void verify_superposition(double **base_phi, double **single_phi, double **super_phi, int n_electrode, int R, int C) { // 初始化叠加数组为零 for (int i = 0; i < R; i++) for (int j = 0; j < C; j++) super_phi[i][j] = 0.0; // 循环处理每个电极 for (int k = 1; k <= n_electrode; k++) { // 设置电极电位：只保留第k个电极的电压，其余电极置零 set_electrode_potential(k); // 调用SOR迭代模块 sor_solve(single_phi, R, C, EPS, MAX_ROUND); // 累加到叠加数组 for (int i = 0; i < R; i++) for (int j = 0; j < C; j++) super_phi[i][j] += single_phi[i][j] * V[k]; } // 对比叠加结果与直接计算结果 double max_diff = 0.0, avg_diff = 0.0, rms_diff = 0.0; for (int i = 0; i < R; i++) { for (int j = 0; j < C; j++) { double diff = fabs(base_phi[i][j] - super_phi[i][j]); max_diff = fmax(max_diff, diff); avg_diff += diff; rms_diff += diff * diff; } } avg_diff /= (R * C); rms_diff = sqrt(rms_diff / (R * C)); printf("最大绝对残差: %f, 平均残差: %f, 均方根残差: %f\n", max_diff, avg_diff, rms_diff); }

这段代码里有一个隐含的“坑”：set_electrode_potential(k)函数需要重新初始化电位矩阵，把非当前电极的电压全部清零。如果只是把电极电位值改掉，而没有清零边界条件，之前轮次的残余电位会污染下一次迭代的初值，导致叠加结果完全错误。我一般会在set_electrode_potential里显式执行一次全矩阵初始化，把边界条件重新赋值一遍，确保每次迭代都在同一套边界设置下进行。叠加定理验证最大的意义在于：它能以极低成本检查整个程序的边界处理和迭代逻辑是否自洽。如果叠加验证的最大绝对残差超过预设的迭代精度一个数量级，那几乎可以肯定边界封闭或者迭代更新顺序里有隐藏的bug，这种bug在单一组的迭代结果里很难被肉眼发现。

5.3 误差来源拆解：截断误差、精度误差和网格离散误差

任何一种数值计算方法的误差都不能被完全消除，只能被控制在可接受范围内。这个程序里的误差来源主要有三类。第一类是计算机本身的截断误差，由于CPU位数有限，浮点数在存储时只能保留有限的有效位数，每一次加减乘除都会丢失最末尾几个bit的精度。这类误差会随着迭代次数累积，但累积速度很慢，在精度要求为0.00002V的数量级下完全可以忽略。第二类是给定精度带来的截断误差，迭代结束条件是前后两次残差小于预设的ε，这意味着所有小于ε的电位波动都被视为“已收敛”而忽略，即使它可能是一个真实存在的微弱物理信号。第三类是网格离散误差，这是数值方法的核心误差来源。网格永远是离散的，电场在网格点之间的连续变化被线性近似替代，网格越粗，误差越大。程序用行列扫描加线性插值来画等势线，本身就是对离散误差的一种补救，但插值阶数只有一阶，在曲率大的地方（比如电极附近）仍然会有明显的位置偏移。

6. 二次开发思路：从圆柱形像管到任意旋转对称结构的改造

拿到这个程序后，大部分人第一反应就是“怎么改代码能自定义电极形状”。原始程序假设了圆柱形像管，所有电极都是完整的环形，电极内孔径半径相同，厚度相同。如果要做更复杂的结构和参数，可以从以下三个方向入手。

第一个方向是网格分布优化。把轴向网格改成非均匀分布，在电极边缘和阴极附近加密网格，在其他区域适当放宽网格间距。这需要修改初始化模块中的步长数组计算逻辑，并且让差分系数计算函数支持任意疏密分布的网格。你可以对比一下在相同最大网格数的条件下，均匀网格和非均匀网格算出来的轴上电位曲线谁更平滑。

第二个方向是算法升级。把连续超松弛迭代换成多重网格法。原始报告参考文献里也提到了多重网格法的应用，它在三维像管电场计算中能明显加速收敛。实现思路是：先在最细的网格上迭代几次，然后限制到粗网格上继续迭代，得到修正项后再延拓回细网格。这个过程循环往复，能大幅度消除低频误差分量。相比单纯调ω，多重网格的收敛速度提升是指数级的。

第三个方向是插值算法替换。程序目前用拉格朗日三项二次插值来计算对称轴上任意点的电位，优点是代码简洁、容易实现，缺点是精度受限。你可以把它替换成三次样条插值，在已知网格点电位和电极边界值的前提下，样条插值能保证二阶导数连续，轴上电位曲线会更光滑，计算电子轨迹时得到的偏导数也更可信。

从投入产出比来看，我建议你优先做第一项改造，即非均匀网格，因为它的工程量最小，但对计算精度和速度的提升立竿见影。网格分布优化的代码改造点集中在初始化模块里，只要把等间距的循环改为按比例函数分布即可。从那以后，我每次把物理模型转成代码时都会强制走一遍“网格独立性验证”——用稀疏网格和加密网格各算一遍，如果结果差异在可接受范围内，才敢把程序交给别人用。希望帮到你。

本文还有配套的精品资源，点击获取