第一章 非线性代数方程组的数
值解法
1.1 直接迭代法
1.2 牛顿法和修正牛顿法
1.3 拟牛顿法
1.4 增量方法
1.5 增量弧长法
2000.3 哈尔滨工业大学 王焕定教授制作 2
非线性问题可分为三类,材料非线性
不管那类非线性问题，最终都归结为一组非
线性方程 Ψ(a)=0,a为待求的未知量。
对许多问题，用某些方法可将 Ψ(a)=0改造成
Ψ(a) =P(a)-R=K(a) a -R=0
的形式。
对非线性问题的方程 Ψ(a)=0，一般只能用数
值方法求近似解答。
、几何
非线性 和边界非线性。 我们只讨论前两类问题。
其 实质是，用一系列线性
方程组的解去逼近所讨论非线性方程组的解。
本章将简单介绍有限元分析中常见的各种求解
非线性方程组的数值方法。
2000.3 哈尔滨工业大学 王焕定教授制作 3
1.1 直接迭代法
当用某些方法将 Ψ(a)=0改造成迭代格式
Ψ(a) =P(a)-R=K(a) a -R=0
后
a1= K(a0)-1R
如果问题是收敛的,a1将比 a0有所改善。
an+1= K(an)-1R Δan=an+1- an
当设范数为 in aa ?? m ax?
或设范数为 2/1T ])[( nnn aaa ??? ?
收敛条件则为 10 ??? ??? nn aa
，设一初始未知量 a0，则由它可得
如此反
复迭代可得
2000.3 哈尔滨工业大学 王焕定教授制作 4
如果考虑到每步迭代
Ψ(an) =P(an)-R=K(an) an -R≠0
将 Ψ(an)视为不平衡力（或失衡力）并作为衡量
收敛的标准
应指出的是，对单变量情况，如讲义图示，
直接迭代实质是“割线”法
10 )( ??? ?? Ra Ψ n
1.1 直接迭代法
返首页
，则收敛条件也可改为
，一定条件下这种
迭代过程是收敛的,但对多自由度情况，由于
未知量通过矩阵 K(an)的元素互相耦合，在迭代
过程中往往出现不稳定现象。
2000.3 哈尔滨工业大学 王焕定教授制作 5
1.2 牛顿法和修正牛顿法
如果将非线性方程 Ψ(a) =0在 an 附近展开，则
又如果 [Ψ’(a)]n的逆存在，则 Δan 近似等于
记 KT(an)=[Ψ’(a)]n,Pn =Ψ(an)
Δan≈-[Ψ’(a)]n-1Ψ(an)
则 Δan≈-KT(an)-1 Pn, an+1=an+Δan
切线矩阵 不平衡力
如此逐步计算，即可得到非线性方程的解答，
这就是牛顿 -拉夫森法。
)1,2,,( 0)(,)()( ??????????? jiaaΨaΨaΨ jnjiniji j ?
Ψ(a) =Ψ(an)+ [Ψ’(a)]n Δan+。。。 =0
或用求和约定可写为
2000.3 哈尔滨工业大学 王焕定教授制作 6
1.2 牛顿法和修正牛顿法
牛顿法要每步都计算切线矩阵 KT（也称刚度）
并解线性方程组，虽精度高，但工作量也大。
)())(,( nkjijnnkjii aΨaΨa 1???? ???
其中 μn的作用是改变切线矩阵 KT的主对角元素，
使奇异性或病态得到改善。 更多的改进方法可
参看沈聚敏, 钢筋混凝土有限元与板壳极限分
析, 等。
此外，在某些非线性问题 (如理想塑性和软化
塑性问题 )中用牛顿法，迭代过程中切线矩阵可
能是奇异的或病态的,为了克服这一现象，可
有多种处理方法，其一是按下式来求
2000.3 哈尔滨工业大学 王焕定教授制作 7
1.2 牛顿法和修正牛顿法
如果在计算的每一步内，矩阵 KT都用初始近
似解 KT0计算，在这种情况下，仅第一步迭代需
要完全求解一个线性方程组，如果将 KT0三角分
解并存储起来，而以后各步迭代中采用公式
10T ))(( ??? aKa n? Ψ )( na
则只需对上式右端项中的 进行回代就行
了。这种方法称为修正的牛顿法。
Ψ )( na
为了提高修正牛顿法的收敛速度可采用某些
过量修正技术。讲义上作了简要介绍，请大家
自己看。
返首页
2000.3 哈尔滨工业大学 王焕定教授制作 8
1.3 拟牛顿法
拟牛顿法的主要思想是：
首先设 (KT)n+1可写成如下修正形式
接着设 (KT)n+1必须满足如下所谓拟牛顿方程
)()()( nnnnn aΨaΨaaK ???? ??? 111
(KT)n+1=(KT)n+(Δ KT)n
式中 (Δ KT)n称为修正矩阵。
由此可建立拟牛顿法迭代格式 (略去了下标 T)
ninini aaa ???? 1)()( nkjnijni aΨKa 1????
)()()( nininnn kkjjij aΨaΨaaK ???? ??? 111
nnn KKK ???? 1
2000.3 哈尔滨工业大学 王焕定教授制作 9
要用拟牛顿法，还需给出修正矩阵的计算。
推导修正矩阵算式的思路是：
(un)和 (vn)是秩 1（或秩 2，讲义为秩 2）的列向
量，将修正矩阵 代入 拟牛顿方程可得
设 (Δ KT)n=(un)(vn)T
如果取 (vn)=(Δ a)n，则当 (Δ a)n≠(0)时
[(KT)n +(un)(vn)T](Δ a)n=(ΔΨ )n
(un)=[(vn)T(Δ a)n]-1[(ΔΨ )n-(KT)n(Δ a)n]
假设 (vn)T(Δ a)n≠0，则有
(Δ KT)n =[((Δ a)n)T(Δ a)n]-1[(ΔΨ )n-
(KT)n(Δ a)n] (Δ a n)T
当 (Δ a)n=(0)时，迭代已收敛,(Δ KT)n =(0)。
2000.3 哈尔滨工业大学 王焕定教授制作 10
讲义上的内容比这里说明的多，但基本思路
是一样的。关于秩 2的算法，请大家自己看。
1.设 (a)0求 (KT)0 ；
对秩 1算法来说，实际使用的步骤为：
2.求 (Δ a)0=-[(KT)0]-1(Ψ )0； a1=a0+Δa0
5.计算 (KT)1=(KT)0+(Δ KT)0;(KT)0= (KT)1.
3,计算 (ΔΨ )0； a0=a1 。
6.重复第 2步，直到达到精度要求为止。
4,计算 (Δ KT)0 = [(ΔΨ )0-(KT)0(Δ a)0]
(Δ a0)T/(Δ a0)T(Δ a)0；
讲义上的塞尔曼公式可用逆矩阵定义验证。
2000.3 哈尔滨工业大学 王焕定教授制作 11
返首页
讲义上给出了
三种方法的对比，
指出了选用什么算法
应考虑的因素。
2000.3 哈尔滨工业大学 王焕定教授制作 12
求解非线性方程组的另一类方法是增量方法。
当问题的性质与加载的历史有关时，如弹塑
性问题，则必须采用增量方法。
1.4 增量方法
使用这种方法需要知道, 荷载, 项 (R)为零时问
题的解 (a)0。在实际问题中,(R)经常代表真实荷
载,(a)0 代表结构位移。在问题的初始状态，
它们均为零。 这种从问题的初值开始，随着荷
载列阵 (R)按增量形式逐渐增大，研究 (a)i的变
化规律的方法，称为增量方法。
2000.3 哈尔滨工业大学 王焕定教授制作 13
设, 荷载, (R)在任一增量步的值为 λ (R)，
λ 为荷载增量因子,(R)为标准荷载列阵,则非
线性方程 (Ψ(a))= (0)可写为
引入切线矩阵且略去高阶小量后可改写为
0???? ?? Δ,Δ,),( ?? ?ΨaΨaΨ a
(Ψ(a,λ ))=(P(a))-λ (R)=(0)
若 λ +Δλ 时的解答为 (a)+Δ (a)，象牛顿法一
样，将 (Ψ((a)+Δ (a)),λ +Δλ ) 按 Taylor级数
展开,则可得
Δ),(Δ T ?? RaKa 1??
2000.3 哈尔滨工业大学 王焕定教授制作 14
设荷载增量因子 λ 分别取 如下值
(a)m+1=(a)m+Δ (a)m
则荷载 (R)可分成 M级，第 m级荷载为 λm(R)，其
增量为 (λm+1-λm)(R)=Δλm(R)。
由此可得
Δ (a)m=[KT((a)m,λ m)]-1Δλm(R)
10 210 ????????? M????
但是，这样做的每一步都将产生误差，结果使
解答漂移。讲义上简单介绍了四种解决漂移的
方法，下面仅对混合法作简单说明，其他方法
请大家自行阅读。
2000.3 哈尔滨工业大学 王焕定教授制作 15
所谓混合法是指，在增量法每一增量步进行
自修正的迭代计算。其 m增量 步 n次迭代的计算
公式为
在实际计算中，对于 m＜ M-1的各增量步的
计算，可以只进行少许几次 (例如 3次 )迭代，
而对于 m=M-1，即最后的一个荷载增量，需耍
使用较多次迭代，以使近似解更接近于真解。
返首页
)()( 1T nmmn mnm PRKa ?? ? ???
自修正
nmnmnm aaa ?????11 ),( mnmnm ?aΨP ?
不平衡力
用混合法求解时，所选取的荷载增量的步长
可以比普通增量算法的步长大一些。
2000.3 哈尔滨工业大学 王焕定教授制作 16
用迭代法或增量法进行极限分析时，在极值
点附近往往可能不收敛。这时可用增量弧长法
来解决 。
为便于理解，以杆单向拉伸为例加以说明。
1.5 增量弧长法
增量弧长法的基本思想是：将 λ作为独立变量，
在每个增量步进行自修正法平衡迭代，在迭代
过程中自动控制荷载因子 λ的取值。也即
)()( 1T nmnmn mnm PRKa ?? ? ???
nmnmnmnm uaaaa ?? ????? 1 ),( nmnmnm ?aΨP ?
前步结果
本步 n次增量
2000.3 哈尔滨工业大学 王焕定教授制作 17
如图所示，矢径可表达为
Rimaimi jiur
??? ?? ??
u
u
u
Rimaimi jiur
??? 111 ??? ?? ??
Rimaimiii jiarrr
????? 111 ??? ???? ????
因为
21 lrr i
m
i
m ????
? ???
由于弧长法引入了如下约束方程
mimim aau ?? ? 1?因此
mimim aau ?? ?? 21?
11 ?? ?? imimim auu ???
由此可得
2000.3 哈尔滨工业大学 王焕定教授制作 18
由矢量代数和约束方程可得
222 )()( lurrr i
m
i
m
iii ??? ????? ???
21 )()( lrrrrr iiiii ??? ?????? ?????
aimimRimaim iuajia
??? )2[()( 111 ????? ??? ???
因此
若记
0)2( ??? iii rrr ??? ??也即
0)])2( 1 ??? ? Rimim j????
RK i mi 1T1 )( ??? imi mi PK 1T2 )( ????
121111 ???? ?? iimiima ?????
则
)()( 1T nmnmn mnm PRKa ?? ? ???
2000.3 哈尔滨工业大学 王焕定教授制作 19
02)( 121 ??? ?? cba imim ????
将其代入约束方程，可得
)()(1 11T11 ???? iia ??
式中系数为
])[()( 12T11 imiiim ub ???? ??? ??
]2)[()( 12T12 imii uc ??? ?? ??
上述式子是从简单情况推出的，如果除 外
均理解为矩阵，即为一般情况的弧长法方程。
im?
一元二次方程有两个根，应取 和 间成
锐角的根。据此可建立判别条件，具体推导这
里从略了。
ir? ir??
2000.3 哈尔滨工业大学 王焕定教授制作 20
0?? ii rr ?? ?
这样做迭代的轨迹很接近圆弧，而计算工作量
可减少很多。
从上述公式可见，求 的工作量是很大的，
为此，可令 和 相互垂直，也即
ir? ir??
1?im??
2000.3 哈尔滨工业大学 王焕定教授制作 21
由相互垂直的条件可得
i
m
ii
m
ii
mi
m u
u
???
??
??
?
?? ?
?
?
1
1
T
1
2
T
1
)(
)(
综上所述，弧长法求解步骤为：
R 1
m??
1)选定荷载参考值,和本步荷载因子,解
得,由 求弧长。1
ma? )()()( 11T11212 ????? ml ?
2)修改切线刚度矩阵并三角化。检查对角元，
正定则加载，负定则加负荷载。若矩阵对应行
列式为零，达到极限荷载。
12?i?
4)由 和 求 和 。R i
mP 11?i?
3)与上一步同时，求不平衡力,),( i
mimim ?aΨP ?
RK i mi 11T11 )( ??? ?? imi mi PK 11T12 )( ??? ???
2000.3 哈尔滨工业大学 王焕定教授制作 22
7)求当前荷载水平和位移。
8)检查是否满足精度要求，满足时加下一级荷
载。步满足重复 2） ~7），直至收敛。
为了加速收敛，减少工作量，国内外学者做
了许多工作，大家如果论文工作需要，可自行
查阅。
利用上述所给的步骤，即可编程实现增量弧
长法计算。
1?im??5)计算荷载因子 ( )。 02)( 121 ??? ?? cba imim ????
1?ima?6)计算 ( )。 1
21111 ???? ?? iimiima ?????
值解法
1.1 直接迭代法
1.2 牛顿法和修正牛顿法
1.3 拟牛顿法
1.4 增量方法
1.5 增量弧长法
2000.3 哈尔滨工业大学 王焕定教授制作 2
非线性问题可分为三类,材料非线性
不管那类非线性问题，最终都归结为一组非
线性方程 Ψ(a)=0,a为待求的未知量。
对许多问题，用某些方法可将 Ψ(a)=0改造成
Ψ(a) =P(a)-R=K(a) a -R=0
的形式。
对非线性问题的方程 Ψ(a)=0，一般只能用数
值方法求近似解答。
、几何
非线性 和边界非线性。 我们只讨论前两类问题。
其 实质是，用一系列线性
方程组的解去逼近所讨论非线性方程组的解。
本章将简单介绍有限元分析中常见的各种求解
非线性方程组的数值方法。
2000.3 哈尔滨工业大学 王焕定教授制作 3
1.1 直接迭代法
当用某些方法将 Ψ(a)=0改造成迭代格式
Ψ(a) =P(a)-R=K(a) a -R=0
后
a1= K(a0)-1R
如果问题是收敛的,a1将比 a0有所改善。
an+1= K(an)-1R Δan=an+1- an
当设范数为 in aa ?? m ax?
或设范数为 2/1T ])[( nnn aaa ??? ?
收敛条件则为 10 ??? ??? nn aa
，设一初始未知量 a0，则由它可得
如此反
复迭代可得
2000.3 哈尔滨工业大学 王焕定教授制作 4
如果考虑到每步迭代
Ψ(an) =P(an)-R=K(an) an -R≠0
将 Ψ(an)视为不平衡力（或失衡力）并作为衡量
收敛的标准
应指出的是，对单变量情况，如讲义图示，
直接迭代实质是“割线”法
10 )( ??? ?? Ra Ψ n
1.1 直接迭代法
返首页
，则收敛条件也可改为
，一定条件下这种
迭代过程是收敛的,但对多自由度情况，由于
未知量通过矩阵 K(an)的元素互相耦合，在迭代
过程中往往出现不稳定现象。
2000.3 哈尔滨工业大学 王焕定教授制作 5
1.2 牛顿法和修正牛顿法
如果将非线性方程 Ψ(a) =0在 an 附近展开，则
又如果 [Ψ’(a)]n的逆存在，则 Δan 近似等于
记 KT(an)=[Ψ’(a)]n,Pn =Ψ(an)
Δan≈-[Ψ’(a)]n-1Ψ(an)
则 Δan≈-KT(an)-1 Pn, an+1=an+Δan
切线矩阵 不平衡力
如此逐步计算，即可得到非线性方程的解答，
这就是牛顿 -拉夫森法。
)1,2,,( 0)(,)()( ??????????? jiaaΨaΨaΨ jnjiniji j ?
Ψ(a) =Ψ(an)+ [Ψ’(a)]n Δan+。。。 =0
或用求和约定可写为
2000.3 哈尔滨工业大学 王焕定教授制作 6
1.2 牛顿法和修正牛顿法
牛顿法要每步都计算切线矩阵 KT（也称刚度）
并解线性方程组，虽精度高，但工作量也大。
)())(,( nkjijnnkjii aΨaΨa 1???? ???
其中 μn的作用是改变切线矩阵 KT的主对角元素，
使奇异性或病态得到改善。 更多的改进方法可
参看沈聚敏, 钢筋混凝土有限元与板壳极限分
析, 等。
此外，在某些非线性问题 (如理想塑性和软化
塑性问题 )中用牛顿法，迭代过程中切线矩阵可
能是奇异的或病态的,为了克服这一现象，可
有多种处理方法，其一是按下式来求
2000.3 哈尔滨工业大学 王焕定教授制作 7
1.2 牛顿法和修正牛顿法
如果在计算的每一步内，矩阵 KT都用初始近
似解 KT0计算，在这种情况下，仅第一步迭代需
要完全求解一个线性方程组，如果将 KT0三角分
解并存储起来，而以后各步迭代中采用公式
10T ))(( ??? aKa n? Ψ )( na
则只需对上式右端项中的 进行回代就行
了。这种方法称为修正的牛顿法。
Ψ )( na
为了提高修正牛顿法的收敛速度可采用某些
过量修正技术。讲义上作了简要介绍，请大家
自己看。
返首页
2000.3 哈尔滨工业大学 王焕定教授制作 8
1.3 拟牛顿法
拟牛顿法的主要思想是：
首先设 (KT)n+1可写成如下修正形式
接着设 (KT)n+1必须满足如下所谓拟牛顿方程
)()()( nnnnn aΨaΨaaK ???? ??? 111
(KT)n+1=(KT)n+(Δ KT)n
式中 (Δ KT)n称为修正矩阵。
由此可建立拟牛顿法迭代格式 (略去了下标 T)
ninini aaa ???? 1)()( nkjnijni aΨKa 1????
)()()( nininnn kkjjij aΨaΨaaK ???? ??? 111
nnn KKK ???? 1
2000.3 哈尔滨工业大学 王焕定教授制作 9
要用拟牛顿法，还需给出修正矩阵的计算。
推导修正矩阵算式的思路是：
(un)和 (vn)是秩 1（或秩 2，讲义为秩 2）的列向
量，将修正矩阵 代入 拟牛顿方程可得
设 (Δ KT)n=(un)(vn)T
如果取 (vn)=(Δ a)n，则当 (Δ a)n≠(0)时
[(KT)n +(un)(vn)T](Δ a)n=(ΔΨ )n
(un)=[(vn)T(Δ a)n]-1[(ΔΨ )n-(KT)n(Δ a)n]
假设 (vn)T(Δ a)n≠0，则有
(Δ KT)n =[((Δ a)n)T(Δ a)n]-1[(ΔΨ )n-
(KT)n(Δ a)n] (Δ a n)T
当 (Δ a)n=(0)时，迭代已收敛,(Δ KT)n =(0)。
2000.3 哈尔滨工业大学 王焕定教授制作 10
讲义上的内容比这里说明的多，但基本思路
是一样的。关于秩 2的算法，请大家自己看。
1.设 (a)0求 (KT)0 ；
对秩 1算法来说，实际使用的步骤为：
2.求 (Δ a)0=-[(KT)0]-1(Ψ )0； a1=a0+Δa0
5.计算 (KT)1=(KT)0+(Δ KT)0;(KT)0= (KT)1.
3,计算 (ΔΨ )0； a0=a1 。
6.重复第 2步，直到达到精度要求为止。
4,计算 (Δ KT)0 = [(ΔΨ )0-(KT)0(Δ a)0]
(Δ a0)T/(Δ a0)T(Δ a)0；
讲义上的塞尔曼公式可用逆矩阵定义验证。
2000.3 哈尔滨工业大学 王焕定教授制作 11
返首页
讲义上给出了
三种方法的对比，
指出了选用什么算法
应考虑的因素。
2000.3 哈尔滨工业大学 王焕定教授制作 12
求解非线性方程组的另一类方法是增量方法。
当问题的性质与加载的历史有关时，如弹塑
性问题，则必须采用增量方法。
1.4 增量方法
使用这种方法需要知道, 荷载, 项 (R)为零时问
题的解 (a)0。在实际问题中,(R)经常代表真实荷
载,(a)0 代表结构位移。在问题的初始状态，
它们均为零。 这种从问题的初值开始，随着荷
载列阵 (R)按增量形式逐渐增大，研究 (a)i的变
化规律的方法，称为增量方法。
2000.3 哈尔滨工业大学 王焕定教授制作 13
设, 荷载, (R)在任一增量步的值为 λ (R)，
λ 为荷载增量因子,(R)为标准荷载列阵,则非
线性方程 (Ψ(a))= (0)可写为
引入切线矩阵且略去高阶小量后可改写为
0???? ?? Δ,Δ,),( ?? ?ΨaΨaΨ a
(Ψ(a,λ ))=(P(a))-λ (R)=(0)
若 λ +Δλ 时的解答为 (a)+Δ (a)，象牛顿法一
样，将 (Ψ((a)+Δ (a)),λ +Δλ ) 按 Taylor级数
展开,则可得
Δ),(Δ T ?? RaKa 1??
2000.3 哈尔滨工业大学 王焕定教授制作 14
设荷载增量因子 λ 分别取 如下值
(a)m+1=(a)m+Δ (a)m
则荷载 (R)可分成 M级，第 m级荷载为 λm(R)，其
增量为 (λm+1-λm)(R)=Δλm(R)。
由此可得
Δ (a)m=[KT((a)m,λ m)]-1Δλm(R)
10 210 ????????? M????
但是，这样做的每一步都将产生误差，结果使
解答漂移。讲义上简单介绍了四种解决漂移的
方法，下面仅对混合法作简单说明，其他方法
请大家自行阅读。
2000.3 哈尔滨工业大学 王焕定教授制作 15
所谓混合法是指，在增量法每一增量步进行
自修正的迭代计算。其 m增量 步 n次迭代的计算
公式为
在实际计算中，对于 m＜ M-1的各增量步的
计算，可以只进行少许几次 (例如 3次 )迭代，
而对于 m=M-1，即最后的一个荷载增量，需耍
使用较多次迭代，以使近似解更接近于真解。
返首页
)()( 1T nmmn mnm PRKa ?? ? ???
自修正
nmnmnm aaa ?????11 ),( mnmnm ?aΨP ?
不平衡力
用混合法求解时，所选取的荷载增量的步长
可以比普通增量算法的步长大一些。
2000.3 哈尔滨工业大学 王焕定教授制作 16
用迭代法或增量法进行极限分析时，在极值
点附近往往可能不收敛。这时可用增量弧长法
来解决 。
为便于理解，以杆单向拉伸为例加以说明。
1.5 增量弧长法
增量弧长法的基本思想是：将 λ作为独立变量，
在每个增量步进行自修正法平衡迭代，在迭代
过程中自动控制荷载因子 λ的取值。也即
)()( 1T nmnmn mnm PRKa ?? ? ???
nmnmnmnm uaaaa ?? ????? 1 ),( nmnmnm ?aΨP ?
前步结果
本步 n次增量
2000.3 哈尔滨工业大学 王焕定教授制作 17
如图所示，矢径可表达为
Rimaimi jiur
??? ?? ??
u
u
u
Rimaimi jiur
??? 111 ??? ?? ??
Rimaimiii jiarrr
????? 111 ??? ???? ????
因为
21 lrr i
m
i
m ????
? ???
由于弧长法引入了如下约束方程
mimim aau ?? ? 1?因此
mimim aau ?? ?? 21?
11 ?? ?? imimim auu ???
由此可得
2000.3 哈尔滨工业大学 王焕定教授制作 18
由矢量代数和约束方程可得
222 )()( lurrr i
m
i
m
iii ??? ????? ???
21 )()( lrrrrr iiiii ??? ?????? ?????
aimimRimaim iuajia
??? )2[()( 111 ????? ??? ???
因此
若记
0)2( ??? iii rrr ??? ??也即
0)])2( 1 ??? ? Rimim j????
RK i mi 1T1 )( ??? imi mi PK 1T2 )( ????
121111 ???? ?? iimiima ?????
则
)()( 1T nmnmn mnm PRKa ?? ? ???
2000.3 哈尔滨工业大学 王焕定教授制作 19
02)( 121 ??? ?? cba imim ????
将其代入约束方程，可得
)()(1 11T11 ???? iia ??
式中系数为
])[()( 12T11 imiiim ub ???? ??? ??
]2)[()( 12T12 imii uc ??? ?? ??
上述式子是从简单情况推出的，如果除 外
均理解为矩阵，即为一般情况的弧长法方程。
im?
一元二次方程有两个根，应取 和 间成
锐角的根。据此可建立判别条件，具体推导这
里从略了。
ir? ir??
2000.3 哈尔滨工业大学 王焕定教授制作 20
0?? ii rr ?? ?
这样做迭代的轨迹很接近圆弧，而计算工作量
可减少很多。
从上述公式可见，求 的工作量是很大的，
为此，可令 和 相互垂直，也即
ir? ir??
1?im??
2000.3 哈尔滨工业大学 王焕定教授制作 21
由相互垂直的条件可得
i
m
ii
m
ii
mi
m u
u
???
??
??
?
?? ?
?
?
1
1
T
1
2
T
1
)(
)(
综上所述，弧长法求解步骤为：
R 1
m??
1)选定荷载参考值,和本步荷载因子,解
得,由 求弧长。1
ma? )()()( 11T11212 ????? ml ?
2)修改切线刚度矩阵并三角化。检查对角元，
正定则加载，负定则加负荷载。若矩阵对应行
列式为零，达到极限荷载。
12?i?
4)由 和 求 和 。R i
mP 11?i?
3)与上一步同时，求不平衡力,),( i
mimim ?aΨP ?
RK i mi 11T11 )( ??? ?? imi mi PK 11T12 )( ??? ???
2000.3 哈尔滨工业大学 王焕定教授制作 22
7)求当前荷载水平和位移。
8)检查是否满足精度要求，满足时加下一级荷
载。步满足重复 2） ~7），直至收敛。
为了加速收敛，减少工作量，国内外学者做
了许多工作，大家如果论文工作需要，可自行
查阅。
利用上述所给的步骤，即可编程实现增量弧
长法计算。
1?im??5)计算荷载因子 ( )。 02)( 121 ??? ?? cba imim ????
1?ima?6)计算 ( )。 1
21111 ???? ?? iimiima ?????
