第二章非线性代数方程组的数值解法

第二章非线性代数方程组的数值解法第一页，共二十八页，编辑于2023年，星期四非线性问题可分为三类：材料非线性、几何非线性和边界非线性。我们只讨论前两类问题。不管那类非线性问题，最终都归结为一组非线性方程Ψ(a)=0，a为待求的未知量。Ψ(a)=0可写成平衡方程的形式Ψ(a)=P(a)-R=K(a)a

-R=0对非线性方程Ψ(a)=0，一般只能用数值方法求近似解答。其实质是，用一系列线性方程组的解去逼近所讨论非线性方程组的解。

分段线性法第二页，共二十八页，编辑于2023年，星期四1.1直接迭代法1.2牛顿法和修正牛顿法1.3增量方法1.4增量弧长法第三页，共二十八页，编辑于2023年，星期四1.1直接迭代法Ψ(a)=P(a)-R=K(a)a

-R=0设初始未知量为a0，根据上式有a1=K(a0)-1R如果问题是收敛的，a1将比a0有所改善。如此反复迭代可得an+1=K(an)-1RΔan=an+1-

an当设范数为或设范数为收敛条件则为将Ψ(an)视为不平衡力并作为衡量收敛的标准第四页，共二十八页，编辑于2023年，星期四对于单变量问题的非线性方程，直接迭代法的计算过程如图1所示。图上给出的是和

之间的关系，而不是

和，之间的关系

对单变量情况，直接迭代实质是“割线”法，一定条件下这种迭代过程是收敛的，但对多自由度情况，由于未知量通过矩阵K(an)的元素互相耦合，在迭代过程中可能会出现不稳定现象。a1=K(a0)-1R直至Δan=an+1-

an

满足收敛条件第五页，共二十八页，编辑于2023年，星期四收敛性第六页，共二十八页，编辑于2023年，星期四1.2牛顿法和修正牛顿法如果将非线性方程Ψ(a)=0在an附近展开，则记KT(an)=[Ψ’(a)]n，Pn

=Ψ(an)Δan≈-[Ψ’(a)]n-1Ψ(an)切线矩阵不平衡力如此逐步计算，即可得到非线性方程的解答，这就是牛顿-拉夫森法。Ψ(a)=Ψ(an)+[Ψ’(a)]nΔan+。。。=0又如果[Ψ’(a)]n的逆存在，则Δan近似等于则Δan≈-KT(an)-1Pn，an+1=an+Δan第七页，共二十八页，编辑于2023年，星期四Δan≈-KT(an)-1Pn，an+1=an+Δan直至Δan满足收敛性第八页，共二十八页，编辑于2023年，星期四1.2牛顿法和修正牛顿法如果在迭代计算的每一步内，矩阵KT都用初始近似解KT0计算，在这种情况下，仅第一步迭代需要完全求解一个线性方程组，如果将KT0三角分解并存储起来，而以后各步迭代中采用迭代公式则只需对上式右端项中的进行回代就行了。这种方法称为修正的牛顿法。返首页第九页，共二十八页，编辑于2023年，星期四使用修正的牛顿法求解非线性方程组，虽然每一步迭代所花费的计算时间减少了，但迭代过程的收敛速度也降低了。为了提高修正牛顿法的收敛速度可采用某些过量修正技术。第十页，共二十八页，编辑于2023年，星期四加速技术：修正牛顿法的过量修正技术－搜索办法返首页在算出

后，新的近似解由下式给出

（i=1,2,···,N）其中

是大于1的正数，它称为过量修正因子。确定的一维搜索办法。将看做N维空间中的搜索方向，我们希望在该方向上找到一个更好的近似值，即找到一个式中的最好的值。虽然沿这一方向，不能期望求得精确解，但我们可以迭择因子(在搜索问题中称为步长因子)，使在搜索方向上的分量为零，即上式是一个关于的单变量非线性方程。通常可用一些比较简单的方法来估算出的大小。第十一页，共二十八页，编辑于2023年，星期四加速技术：Aitken加速技术返首页在算出

后，新的近似解由下式给出其中

是大于1的正数，它称为加速因子。第十二页，共二十八页，编辑于2023年，星期四第十三页，共二十八页，编辑于2023年，星期四第十四页，共二十八页，编辑于2023年，星期四1.2牛顿法和修正牛顿法其中μn的作用是改变切线矩阵KT的主对角元素，使奇异性或病态得到改善。。此外，在某些非线性问题(如理想塑性和软化塑性问题、塑性卸载)中用牛顿法，迭代过程中切线矩阵可能是奇异的或病态的，为了克服这一现象，可有多种处理方法，其一是按下式来求第十五页，共二十八页，编辑于2023年，星期四使用某种算法的计算效率，除了与收敛速度有关外，还与每一步迭代所花费的计算量有关。关于每步的计算量，牛顿法最大，而修正牛顿法最小。因此在实际问题的计算中判断使用哪种方法效率较高，往往需要进行数值实验。总的看来，不同的算法可能适用于不同的问题。选用哪种算法，与所研究问题的性质(例如，对线性的偏离程度)、规模(离散的自由度总数)以及容许误差等因素有关。第十六页，共二十八页，编辑于2023年，星期四求解非线性方程组的另一类方法是增量方法。使用这种方法需要知道“荷载”项(R)为零时问题的解(a)0。在实际问题中，(R)经常代表真实荷载，(a)0代表结构位移。在问题的初始状态，它们均为零。这种从问题的初值开始，随着荷载列阵(R)按增量形式逐渐增大，研究(a)i的变化规律的方法，称为增量方法。当问题的性质与加载的历史有关时，如弹塑性问题，则必须采用增量方法。1.3增量方法第十七页，共二十八页，编辑于2023年，星期四设“荷载”(R)在任一增量步的值为λ(R)，λ为荷载增量因子，(R)为标准荷载列阵,则非线性方程Ψ(a)=0可写为引入切线矩阵且略去高阶小量后可改写为Ψ(a,λ)=P(a)-λR=0若λ+Δλ时的解答为a+Δa，象牛顿法一样，将Ψ(a+Δ(a)),λ+Δλ)按Taylor级数展开，则可得第十八页，共二十八页，编辑于2023年，星期四设荷载增量因子λ分别取如下值am+1=am+Δam则荷载(R)可分成M级，第m级荷载为λmR，其增量为(λm+1-λm)R=ΔλmR。由此可得Δam=[KT(am,λm)]-1ΔλmR但是，这样做的每一步都将产生误差，结果使解答漂移。第十九页，共二十八页，编辑于2023年，星期四自修正算法，在增量法每一增量步进行自修正的迭代计算。其m增量步n次迭代的计算公式为返首页自修正不平衡力使用这种改进的算法，对于每一增量步都相当于做一次修正。第二十页，共二十八页，编辑于2023年，星期四所谓混合法是指，在增量法每一增量步进行自修正的迭代计算。其m增量步n次迭代的计算公式为在实际计算中，对于m＜M-1的各增量步的计算，可以只进行少许几次(例如3次)迭代，而对于m=M-1，即最后的一个荷载增量，需耍使用较多次迭代，以使近似解更接近于真解。返首页自修正不平衡力用混合法求解时，所选取的荷载增量的步长可以比普通增量算法的步长大一些。第二十一页，共二十八页，编辑于2023年，星期四用迭代法或增量法进行极限分析时，在极值点附近往往可能不收敛。这时可用增量弧长法来解决。为便于理解，以杆单向拉伸为例加以说明。1.4增量弧长法增量弧长法的基本思想是：将λ作为独立变量，在每个增量步进行自修正法平衡迭代，在迭代过程中自动控制荷载因子λ的取值。也即前步结果本步n次增量第二十二页，共二十八页，编辑于2023年，星期四如图所示，矢径可表达为uuu有

由于弧长法引入了如下约束方程由此可得第二十三页，共二十八页，编辑于2023年，星期四由矢量代数和约束方程可得也即若记因此要求交叉项为零则将其代入约束方程，可得第二十四页，共二十八页，编辑于2023年，星期四式中系数为上述式子是从简单情况推出的，如果除外均理解为矩阵，即为一般情况的弧长法方程。第二十五页，共二十八页，编辑于2023年，星期四这样做迭代的轨迹很接近圆弧，而计算工作量可减少很多。从上述公式可见，求的工作量是很大的，为此，可令和相互垂直，也即第二十六页，共二十八页，编辑于2023年，星期四由相互垂直的条件可得综上所述，弧长法求解步骤为：1)选定荷载参考值，和本步荷载因子，解得，由求弧长。2)修改切线刚度矩阵并三角化。检查对角元，正定则加载，负定则加负