SOR法(逐次超松弛法)模拟器 返回
数值分析

SOR法(逐次超松弛法)模拟器

一维拉普拉斯方程的迭代求解工具。改变松弛系数ω,就能看到它如何通过使高斯-赛德尔法的更新"过度"来加速收敛,以及最优ω与残差衰减的实时效果。

参数设置
一维网格点数 n
内部未知数的个数(不含边界)
松弛系数 ω
ω=1时为高斯-赛德尔法。2以上则发散
左端边界值 φL
右端边界值 φR
最大迭代次数
达到此次数仍未收敛则停止
计算结果
收敛所需迭代次数
松弛系数 ω
最优 ω(理论)
高斯-赛德尔加速倍数
最终残差
收敛/发散
松弛扫描的动画展示

从平坦的初始状态开始,每次SOR扫描时φ的分布都会向直线形的解松弛。当ω>1(超松弛)时,分布会暂时超过目标,随后稳定下来。

残差收敛历史(对数轴)
迭代次数 vs 松弛系数 ω
理论与主要公式

$$\phi_i^{(k+1)}=(1-\omega)\,\phi_i^{(k)}+\omega\cdot\frac{\phi_{i-1}^{(k+1)}+\phi_{i+1}^{(k)}}{2}$$

SOR法的更新公式。每次扫描中,使用最新的相邻值计算高斯-赛德尔值 (φ_{i−1}+φ_{i+1})/2,再用系数 ω 放大后覆盖。当ω=1时,等同于高斯-赛德尔法。

$$\omega_{opt}=\frac{2}{1+\sin\!\big(\pi/(n+1)\big)}$$

一维拉普拉斯问题(内部点 n 个)最优松弛系数的近似公式。n越大,ω_opt越接近 2。为保证稳定性,ω必须小于 2。

$$\frac{d^2\phi}{dx^2}=0 \;\Rightarrow\; \phi_i=\frac{\phi_{i-1}+\phi_{i+1}}{2}, \qquad \phi(x)=\phi_L+(\phi_R-\phi_L)\frac{x}{L}$$

求解的一维拉普拉斯方程、其离散方程及精确解。精确解是连接左右边界值的直线。

SOR法(逐次超松弛法)简介

🙋
我听说CAE分析最后都要求解一个巨大的联立方程组。SOR法是其中的求解方法之一吗?
🎓
完全正确。不管是结构、热、还是流体,偏微分方程在网格上离散后,最终都变成一个未知数成千上万甚至数亿的线性方程组 Ax=b。如果直接用消去法求解,计算量和内存都会爆炸,所以采用"迭代求解法"逐步逼近正解。SOR法是迭代求解法中最经典的代表,也是理解"松弛"概念最直观的方法。
🙋
"松弛"听起来像放松肌肉一样。是说逐步逼近正解吗?
🎓
直觉上差不多。这个工具求解的是一维拉普拉斯方程 d²φ/dx²=0,答案就是连接左右边界值的"一条直线"。每个网格点满足"自己的值=两侧邻点的平均值"。我们一轮轮扫过所有网格点,用两侧邻点的平均值更新每一点。这就是高斯-赛德尔法。经过多轮扫描,最初平坦的分布会慢慢"松弛"成直线。你把左边的ω改成1.0,看看上面的动画。
🙋
试过了。虽然最后确实变成直线,但需要很多轮啊…… SOR的"超松弛"是为了加快这个过程吗?
🎓
就是这样。高斯-赛德尔每步只移动一点点。与其这么保守,不如"既然要往这个方向走,就大胆点"——把更新幅度乘以ω倍来增幅。公式就是 φ_i ←(1−ω)φ_i + ω×(两侧平均)。把ω改成1.5或1.8,看右上的"迭代次数"卡片。会是高斯-赛德尔的几分之几。下面"迭代次数vs ω"的图会显示出尖锐的U形谷。
🙋
明白了。谷底是最快的。但是ω接近2时,突然变成"发散"了。为什么呢?
🎓
凡事都有个度。ω越大越"大胆",但ω≥2时就过了。每轮迭代反而会把误差放大,而不是缩小。这样误差越来越大,最后数值爆炸——这就是发散。理论上,SOR法收敛的范围严格限于 0<ω<2,范围内有个最快的点叫"最优ω"。一维情况下,最优ω=2/(1+sin(π/(n+1)))。反过来,ω<1虽然稳定,但比高斯-赛德尔还慢,这叫"欠松弛"。所以我们要在1到2之间,尤其是谷底附近。
🙋
原来最优ω依赖网格点数n呢。那实际工程中每次都能用这个公式吗?
🎓
遗憾的是,这么漂亮的ω_opt公式只在这种一维拉普拉斯这样的"理想"问题上用得了。现实的CAE中矩阵很复杂,材料参数空间变化,可能还有对流项,最优ω通常是未知的。工程师常常要试几个案例来调教ω,或者凭经验从ω=1.8左右开始。这也是为什么学这个工具很重要——你能靠实际操作理解ω和迭代次数的关系,面对真实问题时就能更有感觉。

常见问题

高斯-赛德尔法将每个网格点的值更新为"最新相邻值的平均值"。SOR法(逐次超松弛法)使用松弛系数ω增幅此更新,形式为 φ_i ←(1−ω)φ_i + ω·(高斯-赛德尔值)。当ω=1时,SOR法等同于高斯-赛德尔法;当1<ω<2时,更新"过度",使误差中长波长分量的衰减更快,收敛所需迭代次数大幅减少。本工具另外计算ω=1的高斯-赛德尔迭代次数,并显示SOR法的加速倍数。
对于内部点n个的一维拉普拉斯方程在均匀网格上的求解,使反复次数最小的最优ω可近似为 ω_opt = 2/(1+sin(π/(n+1)))。网格越细(n越大),ω_opt越接近2。例如n=25时,ω_opt≈1.785。实际应用中,许多问题的ω_opt在1.7~1.95范围内。无论ω小于或大于此值,迭代次数都会增加,形成在ω=ω_opt处有尖锐最小值的U形曲线。本工具的"迭代次数vs ω"图表可清晰展示此形状。
SOR法收敛的必要条件是 0<ω<2。当ω≥2时,迭代矩阵的谱半径≥1,误差在每次迭代中被放大而发散。残差不仅不减,反而指数增长,最终导致数值溢出(无穷大)。而ω<1被称为"欠松弛",虽然收敛,但速度比高斯-赛德尔法还慢。稳定且高速的范围是1<ω<2,特别是在ω_opt附近。本工具将ω≥2或迭代上限内未收敛的情况判定为"发散"。
SOR法在结构分析、热传导、电磁场、流体分析(特别是压力泊松方程)等椭圆型偏微分方程离散化后的大规模稀疏矩阵求解中,历来被广泛应用。在现代大规模CAE中,共轭梯度法(CG)、多重网格法、代数多重网格前处理的Krylov法收敛速度更快,但SOR法至今仍作为这些高级方法的前处理器(光滑剂)而发挥重要作用。由于算法简单,易于实现,SOR法也是学习迭代求解法原理(残差衰减、松弛系数效应)的最佳教材。

实际应用

热传导与扩散问题的稳态分析:稳态温度分布在无热源区遵循拉普拉斯方程(有热源时是泊松方程)。电子设备散热器、建筑保温材料、地下温度场等用有限差分法/有限体积法离散后,都变成"每点的值=周围平均值"型的联立方程,可用SOR法求解。一维情况下解是直线,二维三维也是同样思路。

流体分析中的压力泊松方程:非压缩流的SIMPLE法、分步法等在每个时间步都要求解一个压力泊松方程。这往往占总计算时间的大部分,所以需要收敛快的迭代法。古老的代码用SOR法(或其对称版SSOR),现代代码用多重网格或前处理共轭梯度法,但SOR仍常作为"光滑剂"(平滑高频误差)被嵌入其中。

静电场与静磁场分析:电位满足的拉普拉斯/泊松方程出现在电容器、绝缘设计、电机磁路等电磁场分析中。根据边界给定的电位(第一类边界条件)求内部电位分布——这就是本工具的场景,SOR迭代能高效求解。

迭代法教学与算法验证:SOR法代码极短,ω变化时收敛加速与发散现象一目了然,是数值分析入门的经典教材。开发新求解器时,常先用SOR生成"正确解"与"收敛目标",再用它验证高级方法的正确性。

常见误区与注意事项

最常见的误解是"ω越大越快"。ω增大确实能加大每步的更新量,但反复次数并非单调递减。最优ω(本工具中为 ω_opt=2/(1+sin(π/(n+1))))前后,再大或再小都会增加反复次数,形成尖锐的U形。图表可清楚显示。而且ω≥2时误差放大、发散。"ω接近2就最快"的想法是错的,要瞄准谷底。

其次是"最优ω公式对所有问题都通用"的误认。ω_opt=2/(1+sin(π/(n+1))) 是理想的一维(或正方形域二维)拉普拉斯问题的理论值,非常特殊。现实CAE中域形不规则、材料参数空间非均匀、有对流项等,矩阵谱半径无法解析求得,最优ω也就未知。此时必须试探,或从ω=1.7~1.9的经验值出发。千万别把公式值当绝对真理。

最后,"SOR是万能高速求解器"的幻想要破灭。SOR实装简单,小问题效果好,但随着网格细化(n增大),所需反迭代数也增,大规模问题上效率天花板明显。即便用最优ω,也赶不上多重网格法或前处理共轭梯度法。现代大规模CAE很少单独用SOR当主求解器,而是作为多重网格的"光滑剂",用来快速削弱高频误差。要认清SOR是迭代法的起点,不是终点。

使用指南

  1. 在网格点数范围(nGridNum)内设置10~100个点,定义一维拉普拉斯方程的离散化计算域
  2. 用网格范围(nGridRange)以米为单位指定物理区域长度(示例:0.1m~1.0m)
  3. 分段调整松弛系数ω(omegaNum)在0.5~1.9范围,观察各个ω对应的迭代次数和残差变化
  4. 输入左端(phiLNum)和右端(phiRNum)的电位或势能值,单位为伏或帕斯卡
  5. 点击"执行"开始迭代计算,查看收敛迭代次数、最终残差、理论最优ω与实际对比等信息

具体计算示例

在网格点数50、计算域长1.0m、左端势能100V、右端0V的条件下,用ω=1.0(高斯-赛德尔法)计算,约250次迭代达残差1.0×10⁻⁶而收敛。相同条件改用ω=1.85(接近最优值),仅需约45次迭代收敛,加速倍数达5.5倍。网格点数扩大到100时,最优ω接近1.95,高斯-赛德尔法的相对加速倍数超过10倍。

工程实务注意

  1. 最优松弛系数ω随网格点数N变化,ω_optimal ≈ 2/(1+π/N),离散化愈细则ω应愈接近1.9~1.99
  2. 由于ω≥2.0会导致发散,金属热交换器(网格密集)的温度分布计算须设ω≤1.95
  3. 左右边界值差异大时(如100V与0V),残差衰减可能局部缓慢,应从残差分布图检查对称性
  4. 大规模模拟(N>1000)应从纯SOR法转向前处理共轭梯度法实现更好的可扩展性