第二章非线性代数方程组的数值解法精选PPT.ppt

第二章非线性代数方程组的数值解法第1页，本讲稿共28页 非线性问题可分为三类：非线性问题可分为三类：材料非线性、几何非材料非线性、几何非线性和边界非线性。我们只讨论前两类问题。线性和边界非线性。我们只讨论前两类问题。不管那类非线性问题，最终都归结为一组非线性不管那类非线性问题，最终都归结为一组非线性方程方程(a)=0，a为为待求的未知量待求的未知量。(a)=0)=0可写成平衡方程的形式可写成平衡方程的形式可写成平衡方程的形式可写成平衡方程的形式(a a)=P(a a)-R=K(a a)a a-R R=0 对非线性方程对非线性方程(a)=0，一般只能用数值方法求近似，一般只能用数值方法求近似解答。其解答。其实质是，用一系列线性方程组的解去逼近所实质是，用一系列线性方程组的解去逼近所实质是，用一系列线性方程组的解去逼近所实质是，用一系列线性方程组的解去逼近所讨论非线性方程组的解。讨论非线性方程组的解。讨论非线性方程组的解。讨论非线性方程组的解。分段线性法分段线性法第2页，本讲稿共28页1.1 直接迭代法直接迭代法1.2 牛顿法和修正牛顿法牛顿法和修正牛顿法1.3 增量方法增量方法1.4 增量弧长法增量弧长法第3页，本讲稿共28页1.1 直接迭代法直接迭代法(a)=P P(a)-R=K(a a)a -R=0=0设初始未知量为设初始未知量为a0 0，根据上式有，根据上式有a a1 1=K(a0)-1 1R如果问题是收敛的，如果问题是收敛的，a1 1将比将比a0有所改善。如此反有所改善。如此反复迭代可得复迭代可得an+1=K(an)-1R aR an=a=an n+1-an当设范数为当设范数为或设范数为或设范数为收敛条件收敛条件则为则为将将将将(a(a(a(an n n n)视为不平衡力并作为衡量收敛的标准视为不平衡力并作为衡量收敛的标准视为不平衡力并作为衡量收敛的标准视为不平衡力并作为衡量收敛的标准第4页，本讲稿共28页对于单变量问题的非线性方程，直接迭代法的计算过程对于单变量问题的非线性方程，直接迭代法的计算过程如图如图1 1所示。图上给出的是所示。图上给出的是和 之间的关系，而不是之间的关系，而不是 和，之间的关系之间的关系 对单变量情况，直接迭代实质是对单变量情况，直接迭代实质是对单变量情况，直接迭代实质是对单变量情况，直接迭代实质是“割线割线割线割线”法，一定条件下法，一定条件下法，一定条件下法，一定条件下这种迭代过程是收敛的，但对多自由度情况，由于未知量通这种迭代过程是收敛的，但对多自由度情况，由于未知量通这种迭代过程是收敛的，但对多自由度情况，由于未知量通这种迭代过程是收敛的，但对多自由度情况，由于未知量通过矩阵过矩阵过矩阵过矩阵K K(a an n)的元素互相耦合，在迭代过程中可能会出现不稳的元素互相耦合，在迭代过程中可能会出现不稳的元素互相耦合，在迭代过程中可能会出现不稳的元素互相耦合，在迭代过程中可能会出现不稳定现象。定现象。定现象。定现象。a1=K(a0)-1R 直至直至an=an+1+1-a an n 满足收敛条件满足收敛条件满足收敛条件满足收敛条件第5页，本讲稿共28页收敛性收敛性第6页，本讲稿共28页1.2 牛顿法牛顿法和修正牛顿法和修正牛顿法 如果将非线性方程如果将非线性方程(a)=0在在an 附近展开，则附近展开，则 记记 KT(an)=(a)n，Pn=(a an n)an-(a)n-1-1(a an)切线矩阵切线矩阵切线矩阵切线矩阵不平衡力不平衡力不平衡力不平衡力 如此逐步计算，即可得到非线性方程的解答，这如此逐步计算，即可得到非线性方程的解答，这就是就是牛顿牛顿-拉夫森法拉夫森法。(a a)=()=(an)+(a)n an+。=0又如果又如果(a)n的逆存在，则的逆存在，则an 近似等于近似等于则则 an-KT(an)-1 Pn，an+1=an+an第7页，本讲稿共28页 an-KT T(an n)-1 Pn，an+1=an+an 直至直至 an 满足收敛性满足收敛性第8页，本讲稿共28页1.2 牛顿法和牛顿法和修正牛顿法修正牛顿法 如果在迭代计算的每一步内，矩阵如果在迭代计算的每一步内，矩阵KT都用初始近都用初始近似解似解KT0计算，在这种情况下，仅第一步迭代需要完计算，在这种情况下，仅第一步迭代需要完全求解一个线性方程组，如果将全求解一个线性方程组，如果将KT0三角分解并存储三角分解并存储起来，而以后各步迭代中采用迭代公式起来，而以后各步迭代中采用迭代公式则只需对上式右端项中的则只需对上式右端项中的 进行回代就行进行回代就行了。这种方法称为了。这种方法称为修正的牛顿法修正的牛顿法。返首页返首页返首页返首页第9页，本讲稿共28页使用修正的牛顿法使用修正的牛顿法求解非线性方程组，求解非线性方程组，虽然虽然每一步迭代所花费的每一步迭代所花费的计算时间减少了，但计算时间减少了，但迭代过迭代过程的收敛速度也降程的收敛速度也降低了。低了。为了提高修为了提高修为了提高修为了提高修正牛顿法的收敛速正牛顿法的收敛速正牛顿法的收敛速正牛顿法的收敛速度可采用某些度可采用某些度可采用某些度可采用某些过量过量过量过量修正技术修正技术修正技术修正技术。第10页，本讲稿共28页加速技术：修正牛顿法的过量修正技术加速技术：修正牛顿法的过量修正技术搜索办法搜索办法返首页返首页返首页返首页在算出 后，新的近似解由下式给出 （i=1,2,N）其中 是大于1的正数，它称为过过量修正因子量修正因子。确定确定的的一维搜索办法一维搜索办法。将 看做N 维空间中的搜索方向，我们希望在该方向上找到一个更好的近似值，即找到一个 式中的最好的 值。虽然沿这一方向，不能期望求得精确解，但我们可以迭择因子 (在搜索问题中称为步长因子步长因子)，使在搜索方向上 的分量为零，即 上式是一个关于上式是一个关于 的单变量非线性方程。通常可用一些的单变量非线性方程。通常可用一些 比较简单的方法来估算出比较简单的方法来估算出 的大小。的大小。第11页，本讲稿共28页加速技术：加速技术：加速技术：加速技术：AitkenAitken加速技术加速技术加速技术加速技术返首页返首页返首页返首页在算出 后，新的近似解由下式给出其中 是大于是大于1的正数，的正数，它称它称为为加速因子加速因子。第12页，本讲稿共28页第13页，本讲稿共28页第14页，本讲稿共28页1.2 牛顿法和修正牛顿法牛顿法和修正牛顿法其中其中 n的作用是改变切线矩阵的作用是改变切线矩阵KT的主对角元素，使的主对角元素，使的主对角元素，使的主对角元素，使奇异性或病态得到改善。奇异性或病态得到改善。奇异性或病态得到改善。奇异性或病态得到改善。此外，在某些非线性问题此外，在某些非线性问题(如理想塑性和软化塑性如理想塑性和软化塑性问题、塑性卸载问题、塑性卸载)中用牛顿法，迭代过程中中用牛顿法，迭代过程中切线矩阵切线矩阵可能是奇异的或病态可能是奇异的或病态的，为了克服这一现象，可有的，为了克服这一现象，可有多种处理方法，其一是按下式来求多种处理方法，其一是按下式来求第15页，本讲稿共28页 使用某种算法的计算效率，除了与收敛速度有使用某种算法的计算效率，除了与收敛速度有关外，还与每一步迭代所花费的计算量有关。关于关外，还与每一步迭代所花费的计算量有关。关于每步的计算量，牛顿法最大，而修正牛顿法最小。每步的计算量，牛顿法最大，而修正牛顿法最小。因此在实际问题的计算中判断使用哪种方法效率较因此在实际问题的计算中判断使用哪种方法效率较高，往往需要进行数值实验。总的看来，不同的算高，往往需要进行数值实验。总的看来，不同的算法可能适用于不同的问题。选用哪种算法，与所研法可能适用于不同的问题。选用哪种算法，与所研究问题的性质究问题的性质(例如，对线性的偏离程度例如，对线性的偏离程度)、规模、规模 (离散的自由度总数离散的自由度总数)以及容许误差等因素有关以及容许误差等因素有关。第16页，本讲稿共28页 求解非线性方程组的另一类方法是增量方法。求解非线性方程组的另一类方法是增量方法。使用这种方法需要知道使用这种方法需要知道“荷载荷载”项项(R R)为零时问题为零时问题为零时问题为零时问题的解的解的解的解(a)0 0。在实际问题中，。在实际问题中，。在实际问题中，。在实际问题中，(R)经常代表真实荷载，经常代表真实荷载，(a)0 0 代表结构位移。在问题的初始状态，它们均代表结构位移。在问题的初始状态，它们均为零。这种从问题的初值开始，随着荷载列阵为零。这种从问题的初值开始，随着荷载列阵(R R)按增量形式逐渐增大，研究按增量形式逐渐增大，研究(a)i i的变化规律的的变化规律的方法，称为增量方法。方法，称为增量方法。当问题的性质与加载的历史有关时，如弹塑性问当问题的性质与加载的历史有关时，如弹塑性问题，则必须采用增量方法。题，则必须采用增量方法。1.3 增量方法增量方法第17页，本讲稿共28页 设设“荷载荷载”(R R)在任一增量步的值为在任一增量步的值为(R)，为荷载增量因子，为荷载增量因子，(R)为标准荷载列阵为标准荷载列阵,则非线性方则非线性方程程(a)=0可写为可写为 引入切线矩阵且略去高阶小量后可改写为引入切线矩阵且略去高阶小量后可改写为(a,a,)=P P(a)-R=0若若若若+时的解答为时的解答为时的解答为时的解答为a a+a，象牛顿法一样，将，象牛顿法一样，将(a+(a a),+)按按按按TaylorTaylorTaylorTaylor级数展开级数展开级数展开级数展开，则可得，则可得第18页，本讲稿共28页 设荷载增量因子设荷载增量因子设荷载增量因子设荷载增量因子分别取分别取分别取分别取如下值如下值am+1=am+am则荷载则荷载(R)可分成可分成可分成可分成MM级，第级，第级，第级，第mm级荷载为级荷载为级荷载为级荷载为 mR R，其增量，其增量为为(m+1-m)R R=mR。由此可得由此可得由此可得由此可得am=KT(amm,m)-1-1mR但是，这样做的每一步都将产生误差，结果使解但是，这样做的每一步都将产生误差，结果使解答漂移。答漂移。第19页，本讲稿共28页自修正算法自修正算法，在增量法每一增量步进行自修正的迭，在增量法每一增量步进行自修正的迭，在增量法每一增量步进行自修正的迭，在增量法每一增量步进行自修正的迭代计算。其代计算。其代计算。其代计算。其m增量增量步步步步n次迭代的计算公式为次迭代的计算公式为返首页返首页返首页返首页自修正自修正不平衡力不平衡力使用这种改进的算法，对于每一增量步都相当于做使用这种改进的算法，对于每一增量步都相当于做一次修正。一次修正。第20页，本讲稿共28页 所谓混合法是指，在增量法每一增量步进行所谓混合法是指，在增量法每一增量步进行所谓混合法是指，在增量法每一增量步进行所谓混合法是指，在增量法每一增量步进行自修正自修正的的迭代迭代计算。其计算。其mm增量增量增量增量步步n次迭代的计算公式为次迭代的计算公式为 在实际计算中，对于在实际计算中，对于 mM-1-1的各增量步的计的各增量步的计算，可以只进行少许几次算，可以只进行少许几次(例如例如3 3次次)迭代，而对迭代，而对于于m=M-1，即最后的一个荷载增量，需耍使用较，即最后的一个荷载增量，需耍使用较多次迭代，以使近似解更接近于真解。多次迭代，以使近似解更接近于真解。返首页返首页返首页返首页自修正自修正自修正自修正不平衡力不平衡力 用混合法求解时，所选取的荷载增量的步长可用混合法求解时，所选取的荷载增量的步长可以比普通增量算法的步长大一些。以比普通增量算法的步长大一些。第21页，本讲稿共28页 用迭代法或增量法进行极限分析时，在极值点用迭代法或增量法进行极限分析时，在极值点附近往往可能不收敛。这时可用增量弧长法来解附近往往可能不收敛。这时可用增量弧长法来解决决。为便于理解，以杆单向拉伸为例加以说明。为便于理解，以杆单向拉伸为例加以说明。为便于理解，以杆单向拉伸为例加以说明。为便于理解，以杆单向拉伸为例加以说明。1.4 增量弧长法增量弧长法 增量弧长法的基本思想是：将增量弧长法的基本思想是：将增量弧长法的基本思想是：将增量弧长法的基本思想是：将作为独立变量，在作为独立变量，在每个增量步进行自修正法每个增量步进行自修正法平衡迭代平衡迭代，在迭代过程，在迭代过程中中自动控制荷载因子自动控制荷载因子的取值的取值。也即。也即前步结果前步结果本步本步n次增量次增量次增量次增量第22页，本讲稿共28页 如图所示，矢径可表达为如图所示，矢径可表达为如图所示，矢径可表达为如图所示，矢径可表达为uuu有有有有 由于弧长法引入了如下约束方程由于弧长法引入了如下约束方程由此可得由此可得第23页，本讲稿共28页 由矢量代数和约束方程可得由矢量代数和约束方程可得也即也即若记若记因此因此要求交叉项为零要求交叉项为零则则则则将其代入约束方程，可得将其代入约束方程，可得第24页，本讲稿共28页式中系数为式中系数为上述式子是从简单情况推出的，如果除上述式子是从简单情况推出的，如果除 外均外均理解为矩阵，即为一般情况的理解为矩阵，即为一般情况的弧长法方程弧长法方程。第25页，本讲稿共28页这样做迭代的轨迹这样做迭代的轨迹这样做迭代的轨迹这样做迭代的轨迹很接近圆弧很接近圆弧，而计算工作量可减，而计算工作量可减少很多。少很多。从上述公式可见，求从上述公式可见，求 的工作量是很大的，的工作量是很大的，为此，可令为此，可令 和和 相互垂直，也即相互垂直，也即第26页，本讲稿共28页 由相互垂直的条件可得由相互垂直的条件可得 综上所述，弧长法求解步骤为：综上所述，弧长法求解步骤为：1)选定荷载参考值选定荷载参考值 ，和本步荷载因子，和本步荷载因子 ，解，解得得 ，由，由 求弧长。求弧长。2)修改切线刚度矩阵并三角化。检查对角元，正定修改切线刚度矩阵并三角化。检查对角元，正定则加载，负定则加负荷载。若矩阵对应行列式为零，则加载，负定则加负荷载。若矩阵对应行列式为零，达到极限荷载。达到极限荷载。3)与上一步同时，求不平衡力与上一步同时，求不平衡力 。4)由由 和和 求求 和和 。自修正法平衡迭代自修正法平衡迭代自修正法平衡迭代自修正法平衡迭代第27页，本讲稿共28页5)计算荷载因子计算荷载因子 。6)计算计算 。7)7)求当前荷载水平和位移。求当前荷载水平和位移。求当前荷载水平和位移。求当前荷载水平和位移。8)检查是否满足精度要求，满足时加下一级荷载。检查是否满足精度要求，满足时加下一级荷载。步满足重复步满足重复2）7），直至收敛。），直至收敛。第28页，本讲稿共28页