黄金比例分割法确定对称逐次超松弛迭代法的最佳松弛因子

2．SSOR-PCG法改进的迭代格式
SSOR法在每步迭代中，需计算系数矩阵与方向向量的乘积，下面引用了该法的改进迭
代格式，避免了该乘积的计算，比原迭代格式节省计算量8％—50％。
SSOR法的改进迭代格式[3]为：
置初值 δ 0 , g 0 = Kδ 0 − R, y 0 = W −1g 0 , z 0 = −Vy 0 , d 0 = W −T z 0 , k = 0 ，
-0.943297E-01
0.154 0.106 0.132 0.106 0.024
运行时间
表2
二者运算时间比较
Marc计算
1分32秒
本文算法计算 17秒
通过所列的表1和表2观察可以发现，本文算法所求出的位移与大型商业软件MARC相比 较精度误差很小，具有较高的精度，且在保证精度的基础上计算时间也大大缩短，有很高的
j =1
j =i +1
j =1
j =i +1
i −1
n
n
i −1
n
∑ ∑ ∑ ∑ ∑ > ( λ − 1 + ωi )Iii − λ ( ωi aij + ωi aij ) > ( λ − 1 + ωi )Iii − λ ( ωi aij + ωi aij + ωi aij )
j =1
j =i +1

黄金比例分割法确定对称逐次超松弛迭代法的最佳 松弛因子
张德宣，杜成斌，孙立国
河海大学工程力学系，南京 (210098)
E-mail：zdx.1234@
摘 要：本文提出了将黄金分割法确定松弛因子与对称逐次超松弛法的改进迭代格式相结合 的迭代算法。算法应用黄金比例分割法确定最佳松弛因子，成功的将其与运行速度和效率很 高的对称逐次超松弛法的改进迭代格式相结合，并给出了迭代收敛性证明，编写了相应的程 序，对一实际结构进行的算例计算表明，与大型商业软件的计算结果相比较，本文所提算法 具有精度高，收敛快的优点。 关键词：黄金分割；对称逐次超松弛法；松弛因子
-0.147813E-01
y
-0.942368E-01
-0.942186E-01
x
0.431350E-04
0.432271E-04
y
-0.353197E-01
-0.352766E-01
x
0.329006E-03
0.328835E-03
误 差％ 0.194 0.082 0.061 0.343 0.145
表1 方向
x y x y x
本文算法与MARC计算节点位移比较
MARC 解
本文算法解
-0.315654E-03
-0.315041E-03
-0.151406E-03
-0.151282E-03
-0.130349E-01
-0.130269E-01
-0.354377E-01
-0.353946E-01
-0.148028E-01
(1)
( ) R： α = y k ,Vy k
式中 ( ， ) 为内积表示，下同。
如果α ≤ ε ，则停止，否则
-1-

⎧τ ⎪
k
⎪δ k
=
+1
( y k ,Vy k ) /(d k ,2z k = δ k +τkdk,
− Vd k ),
⎪ ⎪
y
k
+1
=
避免了计算Ad，比原迭代格式可省 (ra − 2)n 次乘法运算， ra 为A的各行非零元素的平均个
数。因此， ra = 3 是改进迭代格式节省计算量约8﹪； 2ra >> 6 时，可节省计算量约50﹪，
各迭代步计算量接近CG（共轭梯度）法的计算量。
3．松弛因子的选取
本文采用上述对称逐次超松弛法的改进迭代格式求解大型线性方程组，将黄金分割法确 定松弛因子与之相结合，松弛因子的取值对迭代公式的收敛速度影响很大，它的好坏直接影
0.019 0.213 0.122 0.052
y
10151
x
y
10251
x
y
注：表中节点位置参见图1.
-0.168963E-03 0.131212E-01
-0.355535E-01 0.149818E-01
-0.943525E-01
-0.168703E-03 0.131073E-01
-0.355067E-01 0.149659E-01
的迭代次数少，则选取(dot1, dot2)作为收敛区间，如果是dot3对应的迭代次数少，则选取(dot3, dot4)作为收敛区间。
(3)在所选取的收敛区间里循环进行上述的两个步骤，直到选取出迭代次数最少时所对应 的最佳松弛因子。
通过本文中所采用的算例验证，本文方法确定松弛因子的加速效果是明显的，计算的时 间相对较短。在SSOR法改进迭代格式中采用黄金比例分割的综合算法就可以根据不同的具 体工程问题确定最佳松弛因子，该算法简单明了，而且效果显著，大大节省了选取最佳松弛 因子的时间，提高了大型有限元方程组计算的收敛速度。
yk
+τ k (d k
+W
−1 (z k
− Vd k )),
⎪⎪β ⎨
k
= ( y k+1,Vy k+1 ) /( y k ,Vy k ),
(2)
⎪z k+1 = −Vy k+1 + β k z k ,
⎪⎪d k+1 = W −T z k+1 ,
⎪k = k + 1, ⎪⎪⎩转到R
从迭代格式可见，由于V是对角阵，计算Vy和Vd为2n次乘法运算。改进的迭代格式由于
会延长许多，是很不经济的。 黄金分割比例法是一种经典的算法，本文即采用此种方法来确定最佳松弛因子，其具体
步骤和主要思想如下： (1)利用优选法思想,在(1,2)之间选取四个点，即： dot1=1，dot2=dot4-0.618(dot4-dot1)，dot3=dot1+0.618(dot4-dot1)，dot4=2 。 (2)分别取dot2与 dot3作为松弛因子代入迭代程序，比较出最少的迭代次数，如果dot2对应
i, j =1
j =1
j =i+1
∑ ∑ ∑ n
> ( λ − 1 + ωi )Iii − λ ωi aii
i =1
= ⎜⎜⎝⎛ λ
−1+ ωi kn
−
λ ωi
⎟⎟⎠⎞
n i =1
aii
=
⎜⎜⎝⎛
λ
(1 −
k nωi ) kn
+
(ωi
−1) ⎟⎟⎠⎞
n i =1
aii
由1 < ω
<
2 ，进而 ωi
−1 >
差缩小10−s 所需最小迭代次数为 k ≥ s ln10 ，称 R(B) = − ln ρ(B) 为迭代法收敛速度。 − ln ρ(B)
n
∑ 定义3[6] A = (aij )n×n ，如果A的元素满足 aij ≥ aij (i, j = 1,2L, n) ，称A为对角占优矩 j =1 j≠i
n
∑ 阵，记为 A ∈ D0 ；而当 aij > aij (i, j = 1,2L, n) 时，称A为严格对角占优矩阵，记为 j =1 j≠i
减少迭代次数，这并无规律，完全取决于计算当时的情况，如结构划分单元的多少及计算精
度等因素都有一定的关系。
5．算例
设一悬臂梁长20m，高5m,厚1m, 如图1所示。材料弹性模量E＝2.1E7，泊松比ν =0.167,
不计自重，上表面受1000 N/m 的均布荷载。网格划分为平面四节点单元，长度方向200等分， 高度方向50等分，则节点总数为10251，单元总数为10000。
4．收敛性证明
定义1[4] 设A∈ C n∗n ， λ1, λ2 ,L, λn 为A的特征值，我们称
-2-

ρ
(
A)
=
max i
λi
(3)
为A的谱半径。
定义2[5] 由 0 < ρ (B) < 1 ，且误差向量 ε (k) = x (k) − x∗ 满足 ε k = B k ε (0) ，推出欲使初始U 为严格对角占优矩阵，所以 A = D − L − U ≠ 0 ，由于 I − ωL 也为对
角阵，从而 I − ωL −1 ≠ 0 。
分析式（7）后半部分，并由 λ ≥ 1 可得
-3-

i −1
n
i −1
n
∑ ∑ ∑ ∑ (λ − 1 + ωi )Iii − λ ωi aij − ωi aij ≥ ( λ − 1 + ωi )Iii − λ ωi aij − ωi aij
λ < 1 ，结合定义1可知谱半径 ρ (Lω ) < 1，则该迭代方法收敛。
本文论证出了改进的SSOR迭代法收敛的一个充分条件，而由引理1并结合定义2可知，
迭代次数与收敛速度成反比，即 ρ (B) 越小，迭代收敛速度R(B)越快。而一般来说迭代矩阵
B的谱半径 ρ(B) 对ω 极为敏感，但是究竟是ω 取值大能减少迭代次数，还是ω 取值小可以
图1 悬臂梁网格剖分示意图
表1为本文算法与大型商业软件MARC计算的节点位移比较，表2为两者运算时间的比较 （采用的计算机配置为：Pentium(R) D CPU 2.80 GHz , 1.00 GB 内存）。
-4-

节点 2
101 201
5126 10052
1．引言
线性代数方程组 Ax = b (其中A为非奇异系数矩阵,b为列向量)的迭代解法通常有