化工系统工程-非线性代数方程组迭代解法-数值计算方法课件

1、化工系统工程第三章 数学模型求解方法授课教师：杜可杰化学化工学院 化工系统工程第三章 数学模型求解方法3.4 非线性代数方程组迭代解法*（不要求，成熟软件）3.4.1 直接迭代法（direct substitution, direct iteration）3.4.2 部分迭代法(partial substitution)3.4.3 牛顿（Newton）算法 F(x)=0道理相似，理论相通，计算复杂！3.4 非线性代数方程组迭代解法*（不要求，成熟软件）3牛顿法也可写成改进：Broyden法；修正的牛顿法等等 p77-78Quasi-Newton Methods避免了牛顿每轮迭代需要求导数和求逆

2、的复杂运算牛顿法也可写成改进：Broyden法；修正的牛顿法等等 3.5 迭代加速技术* (p78,不要求)3.5.1 DEM（Dominant Eigenvalue Method）加速法 3.5.2 王德人加速法(略,p79) 思想与策略:步长和方向的选择-考虑相邻两次迭代的误差向量的任意一种范数比;只利用变量及相邻两次误差计算范数信息.OR适用范围:单调收敛的过程其他：多维wegstein法本质上F(x)=0代表了迭代收敛与改进的定量信息！3.5 迭代加速技术* (p78,不要求)3.5.1 D3.6 方程组求解的切割技术简介 ， 背景 化工流程模拟是一项非常复杂的任务，往往涉及到成千上万

3、个非线性方程，而且这些方程涉及的函数定义域狭窄、间断点多、连续性和可微性很差，各变量之间交互作用很强，因此如简单地套用现成的算法往往导致求解效率低下，甚至根本无法求解。 在化工流程模拟实践中，总是需要先对数学模型进行一定的处理，将复杂问题分解为一系列较简单的问题，然后再根据具体情况，调用适当的算法进行求解。切割技术就是其中最重要的一种。(基于化工系统工程那部分内容?) 3.6 方程组求解的切割技术简介 ， 背景 案例1， f1f3f2f4x2x4x1x3X3*信息流图 案例1， f1f3f2f4x2x4x1x3X3*信息流图 案例2 ， f1，f2，f3为关于x2，x3，x4的线性方程组，迭代

4、策略？f1f2f3f4X1*信息流图 案例2 ， f1，f2，f3为关于x2，x3，x4的线性Remarks：关于方程组切割技术 (p80，p82)， 1.选取部分变量作为“切割变量”, 则整个迭代过程从大局看变为仅对切割变量的迭代过程，使待求解问题的维数降低。2.运用切割法可将难解的多个变量联立求解的问题分解转化为一系列易解的较少变量联立求解的问题。3.切割法的作用是将原迭代问题转化为另一类迭代问题，其本身并非迭代求根算法。因此，使用了切割法后，仍需使用方程(组)迭代求根算法。(简化)4.若将切割变量计算值直接用于下一轮对切割变量的迭代，则实际上对切割变量是在使用直接迭代法。 Remarks

5、：关于方程组切割技术 (p80，p82)， 1， 5. 使用切割技术，将n个未知数的联立求解问题转化为对少数m个切割变量的求解问题，将复杂的问题分解为一系列简单的问题。关键在于：适当或巧妙地选取切割变量与中间变量，并将原方程组适当地分解为中间变量方程组和切割变量方程组；所选切割变量（组）最宜是那些物理意义明确、易给初值并密切相关的一组变量。 （关于切割技术将在第四章结合化工过程模拟系统（序贯模块法）进行介绍。， 5. 使用切割技术，将n个未知数的联立求解问题转化为对少方程组切割技术， 对切割技术运用得当，可从以下许多方面改善计算效率。容易选取初始值，扩大收敛域。往往选取了切割变量的初值后，中间

6、变量的初值就变得非常容易选取，甚或可直接计算出来。提高数值稳定性。如中间变量方程组的形式适当，比如是线性方程组，就极为利于求解，不易发生误差传递和放大的现象。减少内存占用。由于切割后回避了大规模方程组的联立求解，故通常许多大型矩阵和向量的工作单元都省去了。最易被忽视而最重要的是，提高求解过程的适定性。化工过程数学模型的定义域往往非常狭窄、苛刻，各变量之间取值关系受到物理规律的严格约束。所有变量同时迭代极其容易使得迭代中间结果远超出定义域，违背物理规律，因而导致迭代混乱。正确选取切割变量后就能保证后续的计算在数学模型的定义域内进行，迭代中间结果也能保持物理意义，避免迭代失败。方程组切割技术， 对

7、切割技术运用得当，可从以下许多方面改善计方程组切割技术， 提高了算法的收敛性。经过切割后，复杂问题被分解为若干基本简单问题，有利于选取收敛性好的迭代算法，加快计算速度。并且，使用收敛性高的算法，比较容易设置安全而合理的收敛判据和容差。使用切割技术后，整体上或从外层看是仅对切割变量进行迭代，减少了问题的维数，但在求解中间变量方程组时也常常需要迭代。故此时形成外层迭代中嵌套着内层迭代的情况，似乎比较复杂。其实不然。切割法一般都是针对具有明确物理意义的大型复杂问题，巧妙分解原方程组，使得内层迭代适定性好、求解极其顺利，因而效率很高。如此时采用经典的联立求解，是难于解决上面论述的稳定性、收敛域、适定性

8、、收敛速度等几方面问题的。切割技术是适用于大规模复杂问题的求解技术，绝不是用来求解一般问题甚或习题的。许多方法都有类似的特点，适用于复杂问题的解法与适用于简单问题的解法是截然不同的。假如对极其简单的问题也采用切割技术，反而是生搬硬套，事倍功半。 方程组切割技术， 提高了算法的收敛性。经过切割后，复杂问题方程组切割技术， 巧妙运用切割技术的关键在于适当选取切割变量（组）和对应的切割变量方程（组）。当方程组不太复杂时，变量数和方程数都较少，往往容易直接观察出较佳的切割方案。但对于维数稍多的问题，由于“组合爆炸”的缘故，基本上无法直接观察出合理的切割变量及其对应的方程。即使不考虑问题的物理意义，仅从

9、纯数学的角度，实际上是能够给出选择切割变量及其对应方程的算法的，只不过对此类问题，数学上的合理解往往不是唯一的，也难于从数学角度判定何种方案是最佳的。 方程组分解就是将大型复杂方程组转化为若干简单方程组求解的方法。方程组分解的内容大致包括“分块”和“切割”两部分，目的就是解决切割变量选择及对应的方程选择问题。通俗地讲就是要解决由哪些子方程组去求解哪些未知数以及求解的可行次序。 方程组切割技术， 巧妙运用切割技术的关键在于适当选取3.7 双层法(不要求，闪蒸，精馏） 严格模型 简化模型 3.7 双层法(不要求，闪蒸，精馏） 严格模型 简化模型3.8 常微分方程组初值问题与动态模拟（简介） 引入(

10、概念) 欧拉法 龙格库塔法 3.8 常微分方程组初值问题与动态模拟（简介） 引入(概念案例1:连续流动加热水槽动态模拟 符号说明：F,T 流量(kg/hr)与温度(K)，T0 环境温度(K)，C, 物性，比热(kJ/kg.K)与密度(kg/m3)，K, A 传热系数(kJ/h K m3) 与传热面积(底面积)(m2)，L, V 液位(m)与持液体积(m3)，t 时间(h)建模简化假设：K，C，为常数，理想搅拌（集中参数体系），流体不可压缩，忽略设备热容与传递滞后，焓衡算用显热衡算代替。-液位与温度随时间的动态变化关系3.8.1引入(定义)（p92）案例1:连续流动加热水槽动态模拟 符号说明：F

11、,T 流量连续流动加热水槽动态模型 首先考虑质量衡算： ，即 (1)热量衡算(显热衡算)： 即: (2 )利用(2-12)，代入上式，消去 ,整理得： (3)式(1)与式(3即构成了所需的数学模型（机理模型）。此模型描述了液位与温度随时间的动态变化关系。如给定适当初始条件，即可算出液位与温度随时间的变化曲线(初值问题)。列出数学表达式数学推演连续流动加热水槽动态模型 首先考虑质量衡算： 案例2:物料冷却的数学模型 它是含有自变量t时间)、位置函数T（随时间变化的物料温度）、To（环境温度）、k(降温速率）以及温度的一阶导数 ，是一个常微分方程。案例2:物料冷却的数学模型 它是含有自变量t时间)

12、、位置3.8.1定义（p92）一阶常微分方程： （aX b）(X为纯量;y为纯量或向量)数值解法就是寻求方程(1)、（2）的近似解,即求区间（a，b）上各个分点序列Xn，n=1，2，m的数值解yn化工计算中，只有一些特殊形式的F(x,y)才能计算它们的解析解。大多数情况下只能计算数值解。如果x为时间，则表示动态过程。其求解过程即动态模拟。数学上成为”初值”问题 3.8.1定义（p92）一阶常微分方程：化工计算中，只有一些3.8.2欧拉法（Euler）欧拉:18世纪科学界的代表人物，是那个时代的巨人。他是历来最有才华、最博学的人物之一，也是历史上最多产的一位数学家x0 x1 x2 y(x)y2y

13、1y0斜率向前欧拉法改进欧拉法（预测校正法精度较高）yx基本思想：在求解区间上作等距分割，步长为h，用差商代替导数计算常微分方程。作f(x)在xi处的1阶向前差商：dy/dxy(x(i+1)）-y(xi)/hdy/dx=f(x,y)f(xi,yi） 向后欧拉法y(i+1)=y(i)+hf（x(i+1), y(i+1) （p96）预报值y=f(x,y)3.8.2欧拉法（Euler）欧拉:18世纪科学界的代表人物向前欧拉公式 应用实例例1：假定某物体的温度w因自热而产生的热量可以使物体在每秒钟内以4%的速度增长，同时该物体由于散热可使其温度在每秒种内下降100k，则物体温度随时间变化的微分方程：

14、（t以秒为单位) 分别以初始温x(0)=1500k，y(0)=2500k，z(0)=3500k用欧拉公式预测24秒后的物体温度趋势。向前欧拉公式 应用实例例1：假定某物体的温度w因自热而产生的向前欧拉公式 应用实例解： w0分别以x0=1500，y0=2500，z0=3500代入。计算结果见如表所示。图4-3 三种初始值的温度变化曲线 表4-1Note：表与图的使用向前欧拉公式 应用实例解：图4-3 三种初始值的温度变化曲线4.2.1 向前欧拉公式 实例当自热引起物体温度升高的速度小于散热引起温度下降的速度，物体的温度随时间而逐渐降低；当自热引起物体温度升高的速度与散热引起温度下降的速度平衡时

15、，物体的温度保持不变；当自热引起物体温度升高的速度大于散热引起温度下降的速度，物体的温度随时间而升高。在图中L1，L2，L3分别表示初始值3500k，2500k和1500k的三条温度变化趋势曲线。 4.2.1 向前欧拉公式 实例当自热Remarks方 法显式欧拉隐式欧拉改进欧拉法简单精度低稳定性最好精度低, 计算量大精度提高计算量大Remarks方 法显式欧拉隐式欧拉改进欧拉法简单精度低3.8.3龙格库塔法(补充内容) 思想：由改进的欧拉法可以得到一个启示：即选择不同点的近似斜率的近似值的平均值。效果会更好。龙格库塔法把欧拉法进一步发展了，下面以四阶龙格库塔法为例。方法1：求解一阶常微分方程的

16、数值解3.8.3龙格库塔法(补充内容) 思想：由改进的欧拉法可基本原理：龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。对于一阶精度的欧拉公式有： yi+1=yi+h*K1 K1=f(xi,yi) 当用点xi处的斜率近似值K1与右端点xi+1处的斜率K2的算术平均值作为平均斜率K*的近似值，那么就会得到二阶精度的改进欧拉公式： yi+1=yi+h*( K1+ K2)/2 K1=f(xi,yi) K2=f(xi+h,yi+h*K1) 基本原理：龙格-库塔(Rung

17、e-Kutta)方法是一龙格库塔法算法如果在区间xi,xi+1内多预估几个点上的斜率值K1、K2、Km，并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格库塔公式，也就是在工程中应用广泛的经典龙格库塔算法： yi+1=yi+h*( K1+ 2*K2 +2*K3+ K4)/6 K1=f(xi,yi) K2=f(xi+h/2,yi+h*K1/2) K3=f(xi+h/2,yi+h*K2/2) K4=f(xi+h,yi+h*K3) 龙格库塔法算法如果在区间xi,xi+1内多预估几个点上的关于：龙格库塔法龙格-库塔法具有精度高，收敛

18、，稳定（在一定条件下），计算过程中可以改变步长，不需要计算高阶导数等优点，但仍需计算 在一些点上的值，如四阶龙格-库塔法每计算一步需要计算四次 的值，这给实际计算带来一定的复杂性。关于：龙格库塔法四阶精度的计算框图。四阶精度的计算框图。龙格库塔法小程序 #include #include #define f(x,y) (-1*(x)*(y)*(y) void main(void) double a,b,x0,y0,k1,k2,k3,k4,h; int n,i; printf(input a,b,x0,y0,n:); scanf(%lf%lf%lf%lf%d,&a,&b,&x0,&y0,&n);

19、 printf(x0ty0tk1tk2tk3tk4n); for(h=(b-a)/n,i=0;i!=n;i+) k1=f(x0,y0); k2=f(x0+h/2,y0+k1*h/2); k3=f(x0+h/2,y0+k2*h/2); k4=f(x0+h,y0+h*k3); printf(%lft%lft,x0,y0); printf(%lft%lft,k1,k2); printf(%lft%lfn,k3,k4); y0+=h*(k1+2*k2+2*k3+k4)/6; x0+=h; printf(xn=%lftyn=%lfn,x0,y0); 龙格库塔法小程序 运行结果 input a,b,x0,

20、y0,n:0 5 0 2 20 x0 y0 k1 k2 k3 k4 0.000000 2.000000 -0.000000 -0.500000 -0.469238 -0.886131 0.250000 1.882308 -0.885771 -1.176945 -1.129082 -1.280060 0.500000 1.599896 -1.279834 -1.295851 -1.292250 -1.222728 0.750000 1.279948 -1.228700 -1.110102 -1.139515 -0.990162 1.000000 1.000027 -1.000054 -0.861

21、368 -0.895837 -0.752852 1.250000 0.780556 -0.761584 -0.645858 -0.673410 -0.562189 1.500000 0.615459 -0.568185 -0.481668 -0.500993 -0.420537 1.750000 0.492374 -0.424257 -0.361915 -0.374868 -0.317855 2.000000 0.400054 -0.320087 -0.275466 -0.284067 -0.243598 2.250000 0.329940 -0.244935 -0.212786 -0.218

22、538 -0.189482 2.500000 0.275895 -0.190295 -0.166841 -0.170744 -0.149563 2.750000 0.233602 -0.150068 -0.132704 -0.135399 -0.119703 3.000000 0.200020 -0.120024 -0.106973 -0.108868 -0.097048 3.250000 0.172989 -0.097256 -0.087300 -0.088657 -0.079618 3.500000 0.150956 -0.079757 -0.072054 -0.073042 -0.066030 3.750000 0.132790 -0.066124 -0.060087 -0.060818 -0.055305 4.000000 0.117655 -0.055371 -0.050580 -0.051129 -0.0