弧长法算法.pdf

弧长法算法.pdf

《弧长法算法.pdf》由会员分享，可在线阅读，更多相关《弧长法算法.pdf（11页珍藏版）》请在冰点文库上搜索。

高等数值分析课程研究基于预处理技术和弧长法的非线性方程通用求解子程序总结报告基于预处理技术和弧长法的非线性方程通用求解子程序总结报告清华大学土木工程系土博零陆新征指导教师：

江见鲸2000年12月7日一、课题研究背景一、课题研究背景在结构有限元分析中，经常需要求解非线性方程组fxF=）（

（1）这里向量x为有限元各个节点的位移，f为节点外力。

如果材料的力学性能复杂，特别是当有些材料的力学性能与其应力历史有关时，必须求解一系列的方程组iifxF=）（，并要求x和f都不能过大，以达到“跟踪”结构荷载位移变化全过程的要求。

由于实际计算中的无法预知计算结果，所以当计算接近极值点时，往往因为对节点外力f的设定不合适而导致求解失败或结果“跳跃”，即x或f过大，造成精度的丧失乃至结果的错误。

目前结构分析中解决这一问题的方法主要有以下三种形式：

1、不考虑下降段当我们只关心峰值点时，当发现计算接近峰值点后，不断减小加载步长，这样，可以得到一个较精确的峰值点，但无法得到下降段。

2、迭代强制收敛方法。

如果主要想了解的是下降段而不是峰值点的值，则当计算非常靠近极值点时，迭代如果不收敛或收敛很慢。

这时就让它强制收敛，即当迭代次数大于预先设定的最大迭代次数Ncr时,则认为Ncr次迭代的结果就是非线形方程组的解，进行下一轮计算。

这样可能会在峰值点附近精度不高，特别是当加载步长取得比较大的时候。

3、线形搜索法（本方法是从结构分析软件SAP-StructureAnalysisProgram的说明中看到的，具体方法我也不是很了解）当计算接近峰值点时，对每次计算的结果进行一个小的折减，使程序自动适应下降段。

但如果曲线变化过于突然，或加载步长过大，则本方法失效。

高等数值分析课程研究为了达到对结构荷载位移变化的全过程自动跟踪，我设计了一个基于弧长法的非线性方程通用求解子程序。

由于当计算接近极值点时，刚度矩阵往往是病态的，为了提高在极值点附近的精度，在迭代求解过程中，使用了预处理技术，以提高求解的稳定性和精确性。

二、基本原理及方法二、基本原理及方法2.1、利用弧长法实现对求解过程的自动控制。

、利用弧长法实现对求解过程的自动控制。

2.1.1、基本原理、基本原理结构分析的一般方程为iifxF=）（，为了达到自动控制的目的，我在该方程组中增加一个变量，变原方程组为iifxF=）（

（2）这时，方程组共有n+1个未知量：

nRx，而方程只有n个，因此，必须补充条件。

现在令补充条件为222222Sxf=+（3）则S即为附加的控制变量，我们称之为“弧长”。

这时，方程组变为（）（）11+nnRF，通过引入，从而实现了对计算过程的自动控制。

2.1.2、具体方法、具体方法将方程组

（2）写成牛顿迭代形式，即：

kkkkRfxxK+=）（（4）这里）（kxK即为）（kxF，kx为第k次迭代的位移值，k为第k次迭代的值，kR为节点不平衡力。

令=xfr2，如图2所示，则附加条件（3）可以写成Sr=2（5）即在迭代中，i，i+1点沿一圆弧进行，故称为弧长法。

现在的问题是如何确定）1（+i与）1（+iu。

按图2所示，将迭代向量表示为）1（）（）1（+=iiiurr（6）代入（5），得到高等数值分析课程研究（）（）2）1（）（）1（）（Sururiiii=+（7）又因2）（）（Srrii=（8）所以0）2（）（）1（）1（=+iiiruu（9）又因为（）22）1（）1（）1（）1（）1（fxxuuiiTiii+=+（10）（）22）（）（）（）（）（fxxrriiTiii+=（11）最后得到方程0）2（2

（2）

（2）1

（2）1（）（）1（）1（=+fffxxxiiiiiTi（12）将）1（+ix分为两部分：

IIiIiiixxx）1（）1（）1（）1（+=（13）其中：

fxKIii=+）1（）（（14），）（）1（）（iIIiiRxK=+（15）将（13）代入（12），得到关于）1（+i的二次方程：

（）02）1

（2）1（=+caii（16）系数IiTIixxa）1（）1（1+=，）（）（）1（）1（）（iIIiTIiixxxb+=+）2（）（）1（）1（iIIiTIIixxxc+=+这样，解方程（16），就可以得到）1（+i。

为了简化）1（+i的求解过程，在程序编制过程中，我使用了切平面法，即用垂直于）（ir的向量）1（+iu来代替圆弧，即0）1（）（=+iiur（17）写成矩阵形式为022）1（）（）1（）（=+fxxiiiTi（18）高等数值分析课程研究同样将）1（+ix分为两部分，则可以解得22）（）1（）（）1（）（）1（1ffxxxxiIiiIIiTii+=+（19）与式（13），（16）相比，该式简化了很多。

图2、弧长法迭代过程示意图在实际计算中，除了要知道）（i的数值以外，还需要知道刚度矩阵的正定性，如果刚度矩阵不是正定的，则节点荷载iiff=+1。

2.1.2、弧长法的具体实现方法弧长法的具体实现方法在实际计算中，先进行矩阵正定性判别，利用Lancsoz方法，进行两步Lancsoz过程，得到相应的Ritz值1，n。

因为n2，n为刚度矩阵最小的特征值。

如果02，则认为刚度矩阵为非正定。

弧长法程序模块由三个子程序组成：

1：

subroutineArcLeng（TotalP,TotalU,F,ArcL），功能：

弧长法迭代核心程序2：

subroutineJudge（PD,GK），功能：

判断刚度矩阵的正定性3：

subroutineSolve（GK,Gload,Gdisp），功能：

求解迭代方程组fxx

（1）x（i）x（i+1）r

（1）r（i）r（i+1）u（i+1）

（1）（i）（i+1）（i+1）（i）高等数值分析课程研究子程序ArcLeng的流程图如图3所示：

输入第n步加载时的总荷载向量TotalF，总位移向量TotalU，参考加载向量F，控制弧长ArcLSetGK（GK,TotalU）调用外部子程序SetGK，由总位移向量得到当前系统的切线刚度矩阵GK调用子程序Judge，Judge（PD,GK）判断刚度矩阵的正定性是否加荷方向与加荷方向与参考参考方向相同方向相反Solve（GK,Gload,Gdisp）调用子程序Solve，求解位移向量Gdisp将荷载与位移与参考弧长比较uuFFArcLtT+=2）（，修正荷载大小FF=，uu=GetRF（RF,TotalU+U）SetGK（GK,TotalU+U）得到当前刚度矩阵GK，当前结点不平衡力RJudge（PD,GK）并由现在的刚度矩阵正定性判断修正荷载P的方向PGKui111+=RGKui112+=高等数值分析课程研究+=+1212itituuuu，1211+=iiuuu从而得到第I步迭代结束时的荷载和位移增量PPP+=，uuu+=判断是否达到收敛要求这里以为收敛判据否是子程序结束图3：

子程序ArcLeng的流程图子程序Judge使用Lanczos方法来判断刚度矩阵是否正定。

具体做法是：

进行两步Lanczos过程，得到矩阵22RB，再得到B的特征值1，2。

因为n2，n为刚度矩阵最小的特征值。

如果025时，预处理Lanczos的结果精度对普通高斯消去已经没有什么优势了。

因此，我认为基于Lancsoz方法求解的尝试没有取得预想效果，需要用新的方法代替。

2.2.2、A3预处理技术预处理技术取预处理矩阵M=A3，即取原矩阵的三条主对角线元素为预处理矩阵，对原矩阵进行预处理。

为了比较其效果，这里以Hilbert矩阵为例，比较其效果。

方程bHx=，H为Hilbert矩阵，分别用普通高斯消去法和预处理后的消去法得到结果x，记误差为2xAb，结果如表1所示：

表1、预处理效果比较矩阵的阶数未处理结果的误差预处理结果的误差101.164153218269348E-0101.091393642127514E-010202.607703208923340E-0081.117587089538574E-008301.981854438781738E-0069.313225746154785E-008可见，预处理具有一定的效果。

而且此方法的代价并不大，可以作为一个实用方法使用。

2.2.3、线性方程组求解、线性方程组求解由于计算中实际的刚度矩阵都是稀疏的，因此，我在这里直接使用数学软件包IMSL中的稀疏矩阵求解子程序DLSLXG作为预处理后刚度矩阵的求解子程序。

经比较发现，DLSLXG的解题速度比直接编程计算的速度快一倍以上，当矩阵特别大时，效果更加明显。

三、具体算例三、具体算例为了验证该程序的有效性，我特选择了一个经典的非线性结构分析题目，一个串联双弹簧系统进行了分析。

系统如图4所示：

FF1，X1F2，X2图4：

双弹簧系统此系统在数学上为一二元非线性方程组，未知量为X1,X2，刚度矩阵22RGK。

每根弹簧的“力-位移”关系为75.225.25.0020300+=iiiiiiiixxxxxxFF，如图5所示。

高等数值分析课程研究图5：

弹簧的力-位移关系设第一根弹簧的屈服力20102FF=，在第二根弹簧端部施加荷载，得到的荷载（F）-位移（X2）全曲线如图6所示。

由图6可以看出，本系统在求解时，具有相当的复杂性。

出现了4个极值点，且导数符号也变化了4次。

一般的牛顿法或位移控制法均失败或求解效率很低。

但利用本程序。

可以迅速稳定的得到结果。

迭代次数一般只需2-3次（极值点附近需要6-7次）就可以达到10-7精度，因此，本程序具有较好的使用用途。

为了比较初始参考弧长对计算的影响，我又改变弧长进行了一系列的试算，结果如图7所示。

可以看出，当参考弧长大与一定限度时，其对整个求解过程的控制能力减弱，结果出现“跳跃”。

（如图7中3、4、5组曲线）。

因此，选择合适的初始参考弧长，使计算效率和精度上有一个较好的平衡点，对弧长法有重要意义。

图8比较了各个不同的线性方程组求解程序对结果的影响，应用了预处理技术后的求解子程序在峰值点附近有更好的性能，参考弧长对极值附近的控制更加稳定有效。

（注，为了便于比较，将三条曲线绘于一张图内。

为了清楚起见，第二，第三条曲线分别平移了0.25，0.5。

实际上三条曲线几乎完全相同。

）图9为现在常用结构计算程序的结果与本程序计算结果的对比，取相同的加载步长。

可见，本方法对整个过程计算效果最好。

不但完整的反映出了整个变化过程，精度也相当的高。

00.511.522.533.500.511.522.533.5S/S0F/F0高等数值分析课程研究图6：

荷载，位移全曲线图7：

弧长对计算的影响荷载-位移关系00.511.522.533.544.55012345678位移荷载弧长对计算的影响012345678012345678位移力0.10.20.40.50.8高等数值分析课程研究预处理后在极值点性能的改善图8、各种方程求解方法的比较的结果比较图9计算结果对比另外，还有以下点需要说明：

1、这里列出的其他方法的结果，都是用结构计算软件包得到的。

其当迭代计算次数相同时，这些软件包的计算速度比我的程序要高一些（一般用时为我的程序的80%-90%），我认为，这是这些商品软件开发时进行了细致的优化的结果。

2、如果计算时步长取得很小，线形搜索法也可以得到与本程序相同的结果，但是，其计算工作量几乎要大一倍（在计算本题时，为了得到同等精度的结果，本程序用了普通Gauss消元法的结果Lancsoz算法结果预处理后结果本程序计算结果不考虑下降段计算结果迭代强制收敛计算结果线形搜索法计算结果高等数值分析课程研究150步，线形搜索法用了超过300步，计算时间也要长不少）。

3、实际计算中，在结构刚度变化平缓的阶段，可以适当放大计算步长，也就是说，本程序还有进一步优化的潜力。

四、结论与讨论四、结论与讨论1、弧长法作为求解非线性方程的一种有效方法，一直受到人们的关注和重视。

提出了包括球面弧长法，柱面弧长法，切平面弧长法，常数线性化法等许多方法。

本文利用弧长法中较简便的切平面法，结合结构计算问题特点稍加改进，并使之通用化，得到了较好的结果。

2、Lanczos方法尽管在理论上具有较好的数值稳定性，但是，经过我的实践，我感觉此方法在求解线性方程组上的优势并不是很明显。

还需要更强的预处理技术配合，才能充分发挥其优势。

3、但是，从以上分析也可以看出，由于在结构分析中经常要分析矩阵的正定性，因此，利用Lanczos过程求得矩阵最小特征值的近似值，进而判断矩阵的正定性，仍然使一个十分有效的方法。

4、在采用了A3预处理技术和稀疏矩阵求解子程序后，线性方程组求解的速度和精度都得到了一定的提高，5、我们可以认为，本程序的开发是成功的。

在解决复杂问题时，有着比较明显的优势。

如果能够加入对参考弧长的自动选取，则可以进一步提高本程序的效率。