- 现代化工计算(第二版)
- 徐建良等
- 2689字
- 2025-02-18 00:43:26
任务三 求解一元非线性方程
知识目标
掌握一元三次以上非线性方程求根的原理,熟悉求解一元三次以上非线性方程的方法。
能力目标
能用Newton迭代法手工求解一元三次以上非线性方程的根,能用Excel的单变量求解法求一元三次以上非线性方程的根。
化学化工中的许多问题常常归结为解函数方程f(x)=0。若f(x)是一元线性方程式或一元二次方程式,就可用代数方法求解析解。若f(x)是一元三次或高次的代数方程式,求解析解是不可能的,只能用数值方法求近似解。例如,用于处理真实气体的状态方程Redlich-Kwong(RK)方程(以1mol气体为基准)为
其展开式为
显然,式(1-21)为体积V的一元三次非线性方程。如果要计算在一定温度T和压力p时气体的摩尔体积,则需求解f(V)这样一个一元三次非线性方程。
当然,采用数值计算法能处理上述问题,但计算过程往往较繁杂且需花费大量时间,但若以计算机为工具,利用数值计算法处理此问题,就很容易得到极为准确的数值解。本任务主要介绍处理一元非线性代数方程的常用数值计算方法及借助于计算机的具体处理步骤。
一、一元非线性方程求根方法
(一)逐步扫描法
逐步扫描法是一元非线性方程求根初始近似值常用的方法,若f(x)是n次多项式,对应的方程为n次代数方程,这时方程的根也称为多项式的根。根有实根和复根之分,这里仅介绍实根的求法。用数值方法求方程的根可分为两步,先找出根的某个粗略近似值,又称为“初始近似值”,然后再将初始近似值逐步加工成满足精度要求的结果。
设待解方程为
f(x)=0
在直角坐标系中给出相应于y=f(x)的曲线,如图1-28所示。显然,此曲线与x轴的交点就是方程f(x)=0的根。若函数f(x)在区间(a,b)连续,且f(a)与f(b)异号,则区间(a,b)内必定至少有一个实根。若函数f(x)在区间(a,b)连续并单调(单调上升或单调下降),则在区间(a,b)必定只有一个实根(图1-29)。这时,选定一个步长h,并计算函数f(a)、f(a+h)、f(a)与f(a+h)两函数值乘积,若两函数值乘积大于零,说明该区间内无实根。这时,把区间的终点作为下一次迈步的始点。再计算f(a+2h)的值,如此反复向前迈步,直至两函数值的乘积小于或等于零,即直至相邻两个函数值异号,此时可把两函数乘积小于或等于零的区间的始点作为方程式根的近似值。这个方法叫迈步法或逐步扫描法。

图1-28 y=f(x)曲线

图1-29 单调变化的y=f(x)曲线
显然,步长缩小,精度提高。当精度要求较高时,步长要求很小,计算机循环计算的次数将会很大,所耗机时较大。因此,通常不用该方法求根的精确值,只用于求根的初始近似值。
用迈步法求得非线性代数方程根的初始近似值后,可用以下介绍的二分法、牛顿法和弦截法等寻求根的精确值。
☆思考:何谓逐步扫描法?
(二)二分法
设求解方程f(x)=0通过逐步扫描法已知有根区间为(x1,x2),如图1-30所示。现取x1与x2的中点x0,即

图1-30 二分法求根的精确值
x0=(x1+x2)/2
从而将区间(x1,x2)分为相等的两部分。然后检查f(x0)与f(x1)的符号是否相同,如为同号,根必在x0与x2之间,如图1-30(a)。令
x1=x0
如为异号,则根必在x1与x0之间,如图1-30(b)。令
x2=x0
然后再取新区间(x1,x2)的中点,重复以上步骤,直至x1与x2之间的距离小于某指定值ε为止。取前一区间的中点x0作为方程式根的精确值。这个方法称为“二分法”或“对分法”。只要区间(x1,x2)内有根,用二分法一定能够求出结果,因而它是最保险的一种方法,但其最大的缺点是收敛速度较慢。
需要说明的是,用上述方法求根,只能得到方程式的一个实根。实际上,若方程式有数个实根时,只需增加一个终值B,当求出第一个实根之后,检查有根区间的终点x2是否小于终值B。若小于终值B,则再施行迈步法求下一个实根的近似值,进而求出精确值。如此反复进行下去,直至达到或超过终值B为止。
☆思考:何谓二分法?
(三)牛顿法
牛顿法亦称牛顿-拉福森(Newton-Raphson)法。由于这个方法的计算结果颇佳,而计算过程亦较简单,故被普遍地使用。
如图1-31所示,假设方程f(x)=0有一个实根x*。取一初值x0,过x0作x轴垂线交于曲线f(x)于点P0,过P0点作曲线f(x)的切线并于x轴相交于x1点,显然x1点较x0点更接近于根x*,如果|x1-x0|<ε,则方程根x*=x1,否则按上述同样方法过x1作x轴垂线交于曲线f(x)于点P1,过P1点作曲线f(x)的切线并于x轴相交于x2……,直到|xk+1-xk|<ε为止,方程的根为

图1-31 牛顿迭代法
x*=xk+1
由图1-31可得
则
于是可写出牛顿迭代法的通式
式中 xk,xk+1——第k、k+1次求得的方程近似根。
只要函数y=f(x)的导数f'(x)易得,初值适当,用牛顿迭代法求非线性方程f(x)=0的根既快又正确。但是必须注意,起始值若不在根的附近,牛顿迭代公式不一定收敛。所以,在实际使用中,牛顿迭代法最好与逐步扫描法结合起来,先通过逐步扫描法求出根的近似值,然后再用牛顿迭代公式求其精确值,以发挥牛顿法收敛速度快的优点。
☆思考:何谓牛顿迭代法?
(四)弦截法
虽然牛顿迭代法收敛速度快,但它要求计算f'(x)的值。在科学与工程计算中,常会遇到f'(x)不易计算或算式复杂而不便计算的情况。下面介绍一种既具有牛顿迭代法收敛速度快的优点、又不用求导数f'(x)的弦截法。
弦截法的基本思想与牛顿迭代法相似,即将非线性函数f(x)线性化后求解。两者的差别在于弦截法实现函数线性化的手段采用的是两点间的弦线,而不是某点的切线。
如图1-32所示,若已知非线性函数f(x)的根区间(x0,x1),过x0、x1作垂线交函数f(x)于P0、P1点,连接P0、P1交x轴于x2点。

图1-32 弦截法
若f(x2)=0,则方程的解为x*=x2;若f(x2)f(x1)>0,如图1-32(a),则用x2代替x1;若f(x2)f(x1)<0,如图1-32(b),则用x2代替x0;再用新得到的两个点用以上方法继续迭代,直到相邻两次值满足|xk+1-xk|<ε为止,不难可得弦截法的迭代格式为
式中 xk-1,xk,xk+1——第k-1、k、k+1次求得的方程近似根。
与牛顿法只需给出一个初值不同,弦截法需要给出两个迭代初值x0和x1。如果与逐步扫描法结合起来,则最后搜索的区间的两个端点值常可作为初值x0和x1。
☆思考:何谓弦截法?
弦截法虽比牛顿法收敛速度稍慢,但在每次迭代中只需计算一次函数值,又不必求函数的导数,且对初值x0和x1要求不甚苛刻,因此是工程计算中常用的有效计算方法之一。
二、手工求解一元非线性方程
上面讲解了几种求解一元非线性方程的方法,下面通过实例介绍手工求解一元非线性方程的具体步骤。
【例1-8】 某一常压气相反应体系,当反应达平衡时其中某一组分的平衡分压p(单位为atm)符合以下方程:
4p3-1.640p2+1.640p-0.410=0 (1-24)
试求其分压的数值?
解:以Newton迭代法为例
步骤1:令 y=4p3-1.640p2+1.640p-0.410 (1)
则 y'=12p2-3.280p+1.640 (2)
步骤2:确定p的初值,由题意知为常压反应系统,则其分压可取为p0=0.3atm。
步骤3:将p0=0.3atm代入式(1)、式(2)计算得:y0=0.04240、。
步骤4:将y0、代入式(1-22)计算得:
步骤5:计算p0、p1之间的相对误差ε1
步骤6:判断误差ε1是否符合计算精度(一人为规定的两次计算值之间的相对误差值,手工计算时可取10-3,计算机编程计算可取10-4,甚至更小些)要求,若满足精度要求,计算结束,此次计算值即为解,否则将本次计算值作为初值,重复步骤3~6,直到满足计算精度要求。
本次计算的误差ε1>10-3,需重复计算步骤3~6,具体计算结果如表1-21所示,从表中数据可知,通过三次迭代计算,其迭代精度已达10-6数量级,所以,该组分的平衡分压p为:
p=0.2749atm
表1-21 例1-8手工Newton迭代法计算数据汇总表

☆思考:如何利用Excel求解一元非线性方程?
三、采用Excel求解一元非线性方程
用Excel求解一元非线性方程可分为三种方法:①基于Excel单元格的Newton迭代法;②Excel自带的单变量求解法;③规划求解法。以下仍以例1-8为例,详细说明采用以上三种方法如何求解非线性方程的过程。
(一)基于Excel单元格的Newton法
步骤1:将原方程、一阶导数式的各系数输入单元格,如图1-33中单元格D3、E3、F3、G3、E4、F4、G4。

图1-33 例1-8Excel单元格结合Newton迭代法求解
步骤2:第一次迭代取初值为0.3,如图1-33中单元格D6,迭代精度取1×10-4,如单元格F6。建计算过程表格,如单元格A7~G7,在单元格A8~G8中对应分别输入:“1”、“0.3”、“=$D$3*B8 3+$E$3*B8 2+$F$3*B8+$G$3”、“=$E$4*B8 2+$F$4*B8+$G$4”、“=B8-C8/D8”、“=(E8-B8)/B8”、“=IF(ABS(F8)<$F$6”,“满足精度”,“重新计算”)。
在单元格C8中计算得函数值0.0424,单元格D8中得一阶导数值1.736,单元格E8中输入的是牛顿迭代公式(1-22)计算得到第一次迭代结果值0.275576,单元格F8中计算的是第一次迭代结果与假设的初值之间的相对误差值-8.141×10-2,单元格G8用来判断计算是否可以结果,采用的是Excel自带的IF函数,其中ABS为绝对值函数。当单元格G8中判断需“重新计算”时,需进行第二次迭代计算。
步骤3:在单元格A9中输入“2”,B9中输入“=E8”,即将第一次计算的结果作为第二次计算的初值,选中单元格C8~G8,利用Excel的自动填充功能向下拉,让Excel自动填充单元格C9~G9,由图1-33可见第二次迭代的计算结果。当G9中仍为:“重新计算”时,则需进行第三次迭代计算。
步骤4:由于第三次迭代计算方法与前两次完全一样,可以让Excel自动进行。选中单元格B9~G9,利用Excel自动填充得B10~G10的结果,由图1-33可见,第三次迭代计算结果“0.274901”与第二次计算的结果“0.274902”相对误差已达10-6数量级,满足计算精度要求,计算到此结束,计算结果与手工计算完全一致。
若第三次迭代计算精度仍不能满足要求,则需进行第四、五……,方法与步骤4完全相同,直到计算满足精度为止。
注意:此例中为便于使用Excel单元格的自动填充功能,单元中的输入格式有两种:绝对引用格式和相对引用格式,这两种格式的详细区别请参阅Excel相关章节的内容。
(二)Excel自带的单变量求解法
步骤1:打开Excel,选择某一单元格,如图1-34中的单元格B14,输入:

图1-34 例1-8Excel单变量求解步骤1
“=4*A14 3-1.64*A14 2+1.64*A14-0.41”。
此时由于单元A14中没有任何信息,Excel将其视为“0”,输入完毕确定,在单元格B14中将得到原函数常数项数值,即“-0.41”。
单元格A14相当于原表达式中的分压p,若在此单元格中输入某一数据可使单元格B14中的值为0,则此数据为该函数的一个解,单变量求解就据此原理。
步骤2:选中单元格B14,如图1-34,单击“数据”⇒“假设分析”⇒“单变量求解”,出现“单变量求解”对话框(见图1-35),在各空格中依次填入以下信息。

图1-35 例1-8Excel单变量求解步骤2——“单变量求解”对话框
目标单元格:B14;
目标值:0;
可变单元格:$A$14或A14。
步骤3:在图1-35中按“确定”按钮后⇒图1-36“单变量求解状态”对话框。单击确定按钮,答案出现在单元格A14中(如图1-35所示):0.27490124,很显然,此值与前面计算的结果完全一致。

图1-36 “单变量求解状态”对话框
需要说明的是单元格B14中的值为0时,说明此根为一精确解。但有时单元格B14中是一个很小的数值,Excel可能显示的值为“0”而并非真为0。若要减小误差使解更精确,可采取以下方法操作:单击左上角“Excel按钮”⇒“Excel选项”⇒“公式”⇒出现“Excel选项”卡,如图1-37所示,可将“计算选项”中的“最多迭代次数”和“最大误差”值进行修改,如将最多迭代次数由“100”改为“10000”,此值最大值为32767。最小误差由“0.001”改为“0.0001”等,这样可使计算精度得到提高。

图1-37 重新计算选项卡
注意,此处的最小误差不能改得过小,否则有可能导致超出最大循环次数而无法计算。
☆思考:如何利用ExcelVBA编写一段求解一元三次非线性方程的小程序?
(三)规划求解法
步骤1:打开Excel,选择某一单元格如图1-38中的单元格B20,按式(1-24)输入前三项:“=4*A20 3-1.64*A20 2+1.64*A20”。此时由于单元A20中没有任何信息,Excel将其视为“0”,输入完毕确定,在单元格B20中将显示“0”。

图1-38 例1-8规划求解步骤1
步骤2:选中单元格B20,单击“数据”,选择“规划求解”,调出“规划求解参数”对话框,如图1-39所示。选择“值为( <文字颜色>V</文字颜色> )”,输入方程的常数项值“0.41”。可变单元格输入框中输入:$A$20。单击“求解”按钮⇒图1-40的“规划求解结果”对话框。

图1-39 例1-8“规划求解参数”对话框
步骤3:单击图1-40中的“确定”按钮⇒规划求解结果,如图中单元格A20中的数据“0.27490143”,可见计算结果与前面两种方法完全一致。

图1-40 例1-8“规划求解结果”对话框
注意:在Office软件安装时,规划求解往往不会被安装,可以通过加载实现。步骤是:单击软件左上角“Office按钮”⇒“Excel选项”⇒“加载项”⇒选择“规划求解”加载项⇒“加载宏”对话框⇒在“规划求解”项打勾即可。
四、技能拓展——ExcelVBA自定义函数求解一元三次非线性方程
应用牛顿迭代法计算原理,采用VBA自编迭代函数求解,可使求解过程大为简化,以下仍以例1-8为例说明具体步骤。
步骤1:打开Excel⇒“开发工具”⇒“VisualBasic”编辑器⇒插入⇒模块⇒“添加过程”对话框,如处理的是一元三次非线性方程,则可输入函数名“Newton3”,如图1-41所示。

图1-41 例1-8添加自定义VBA函数示意图
步骤2:选择好相关类型后,按“确定”,出现VBA代码编辑窗口,在编辑窗口编写VBA代码,如图1-42所示。

图1-42 例1-8自编VBA代码窗口示意图
步骤3:在Excel表格中分别输入方程初值、迭代精度及方程中各系数之值,如图1-43中单元格D6、F6、D3、E3、F3、G3。

图1-43 例1-8利用自编VBA函数求解示意图
步骤4:在图1-43单元格A17中输入:=Newton3(D6,F6,D3,E3,F3,G3),按“Enter”即得方程的解为:0.27490124。
显然,采用自编VBA函数求解结果与单变量求解结果完全相同。