1、第二講 非線性方程（組）求解與迭代法,隋志軍 化工學(xué)院軟件應(yīng)用教科組 2006-10,本章知識要點(diǎn),數(shù)值計(jì)算 單個(gè)非線性方程求解 非線性方程組求解 迭代法,MATLAB 求解非線性方程（組）的相關(guān)函數(shù),本章的所要解決的典型問題,在945.36kPa（9.33atm）、300.2K時(shí)，容器中充以2mol氮?dú)?，試求容器體積。 已知此狀態(tài)下氮?dú)獾腜-V-T關(guān)系符合范德華方程，其范德華常數(shù)為a4.17atmL/mol2，b0.0371L/mol,數(shù)學(xué)模型：,非線性方程（組）在化學(xué)工程中的作用,多組分混合溶液的沸點(diǎn)、飽和蒸氣壓計(jì)算 流體在管道中阻力計(jì)算 多組分多平衡級分離操作模擬計(jì)算 平衡常數(shù)法求解化學(xué)

2、平衡問題 定態(tài)操作的全混流反應(yīng)器的操作分析,求解非線性方程的方法,二分法 不動點(diǎn)迭代 威格斯坦法迭代 牛頓法 割線法,迭代法原理,不動點(diǎn)迭代：,迭代法意義示意圖,迭代法的收斂性,迭代法意義示意圖,Wegstein法意義示意圖,牛頓法意義示意圖,割線法意義示意圖,多項(xiàng)式函數(shù),Matlab多項(xiàng)式函數(shù) Polyval多項(xiàng)式的值 Polyder多項(xiàng)式微分 Polyfit多項(xiàng)式擬合 Conv，Deconv多項(xiàng)式乘法和除法 Roots多項(xiàng)式求根,多項(xiàng)式：,Matlab表達(dá)多項(xiàng)式方法：,簡單多項(xiàng)式函數(shù)的使用,a=1 11 55 125; b=1 1;1 1; c=polyval(a,b)求多項(xiàng)式a在b處的值

3、 q=2,3,5; aq=conv(a,q)多項(xiàng)式乘法 d=poly2sym(aq)多項(xiàng)式向量表示為符號多項(xiàng)式 e=deconv(aq,a)多項(xiàng)式除法,多項(xiàng)式求根函數(shù)roots,求解方程：,p=2 -5 6 -1 9; sol=roots(p) 結(jié)果： sol = 1.6024 + 1.2709i 1.6024 - 1.2709i -0.3524 + 0.9755i -0.3524 - 0.9755i,P = 9.33; T = 300.2; n = 2; a = 4.17; b = 0.0371; R = 0.08206 coef=P, -(P*n*b+n*R*T), a*n2 , -a*b

4、*n3; V=roots(coef) 結(jié)果： V =5.0028 0.2429 0.1092,非線性方程求解函數(shù)fzero,調(diào)用格式： x,fval,exitflag,output = fzero(fun,x0,options, p1, p2, .) 此函數(shù)的作用求函數(shù)fun在x0附件的零值點(diǎn)，x0是標(biāo)量。 fval函數(shù)在解x處的值 exitflag程序結(jié)束情況，0，程序收斂于解；0，程序沒有收斂；0，計(jì)算達(dá)到了最大次數(shù) output一個(gè)結(jié)構(gòu)體，提供程序運(yùn)行的信息；output.iterations，迭代次數(shù)； output.functions，函數(shù)fun的計(jì)算次數(shù)； output.algor

5、ithm，使用的算法 options選項(xiàng)，可用optimset函數(shù)設(shè)定選項(xiàng)的新值 fun可以是函數(shù)句柄或匿名函數(shù)。,fzero函數(shù)的使用,1) sinx在3附近的零點(diǎn),2) cosx在1,2范圍內(nèi)的零點(diǎn),3),4),fzero(sin,3),fzero(cos,1,2),fzero(x) x3-2*sin(x),1) 或者： function Cha2demo1 x=fzero(fun,1) function y=fun(x) y=x3-2*sin(x);,fzero(x) x3-2*x-5,1); roots(1 0 -2 -5),fzero函數(shù)初值的選取,P = 9.33; T = 300

6、.2; n = 2; a = 4.17; b = 0.0371; R = 0.08206; V=fzero(V) P*V3-(P*n*b+n*R*T)*V2+ a*n2*V -a*b*n3,0),fzero函數(shù)初值的選取,以t為自變量，取值范圍為-10t10，a,b為參數(shù)，本例取值分別為0.1，0.5,的零點(diǎn),function Cha2demo2 a=0.1;b=0.5;t=-10:0.01:10; Y=sin(t).2.*exp(-a*t)-b*abs(t); clf,plot(t,Y,r);hold on,plot(t,zeros(size(t),k); xlabel(t);ylabel(

7、y(t),hold off zoom on n=input(How many zero points are there?); tt,yy=ginput(n);zoom off for i=1:n t0(i),y(i),exitflag=fzero(t) sin(t)2*exp(-a*t)-b*abs(t),tt(i); end disp(The zero points are:) fprintf(%.4ft,t0) fprintf(n),fzero函數(shù)初值的選取,在945.36kPa（9.33atm）、300.2K時(shí)，容器中充以2mol氮?dú)?，試求容器體積。 已知此狀態(tài)下氮?dú)獾腜-V-T關(guān)系符

8、合范德華方程，其范德華常數(shù)為a4.17atmL/mol2，b0.0371L/mol,function PVT P = 9.33; % atm T = 300.2; % K n = 2; % mol a = 4.17; b = 0.0371; R = 0.08206; V0 = n*R*T/P V,fval = fzero(PVTeq,V0,P,T,n,a,b,R) % - function f = PVTeq(V,P,T,n,a,b,R) f = (P + a*n2/V2) * (V-n*b) - n*R*T;,fzero函數(shù)初值對結(jié)果的影響,function Conv=Cha2demo4 T

9、0 = 450; x0 = 0:0.1:1.0; n=length(x0); for i=1:n x(i)= fzero(NonlinEq,x0(i),T0); end disp(The conversion could be) fprintf(%.4ft,x) fprintf(n) % - function f = NonlinEq(x,T0) T = T0 + 250*x; f = 0.25*(1-x)2*exp(20-10000/T) - x;,結(jié)果： The conversion could be 0.04080.04080.28040.28040.28040.28040.83630.

10、83630.83630.83631.0951,fsolve函數(shù),與fzero函數(shù)只能求解單個(gè)方程的根不同，fsolve函數(shù)可求解非線性方程組的解。其算法采用的是最小二乘法。 調(diào)用格式： x,fval,exitflag,output,jacobian = fsolve(fun,x0,options, p1, p2, .) 輸入輸入變量的意義同fzero函數(shù)。輸出變量中的jacobian為函數(shù)fun在x處的Jacobian矩陣。,fsolve函數(shù)的使用,function Cha2demo6 x0=1 1 1; x=fsolve(fun,x0) function y=fun(x) y(1)=sin(

11、x(1)+x(2)2+log(x(3)-7; y(2)=3*x(1)+2x(2)-x(3)3+1; y(3)=x(1)+x(2)+x(3)-5; 解得結(jié)果如下：x0.5991 2.3959 2.0050,fsolve函數(shù)的應(yīng)用,在銅管內(nèi)在1 atm下將異丙醇加熱到400。已知銅是生產(chǎn)丙酮和丙醛的催化劑，或許還有某些異丙醇異構(gòu)化為正丙醇。這三種產(chǎn)物的生成可用如下三個(gè)獨(dú)立反應(yīng)表示： iC3H7OH(IP) = n C3H7OH(NP) K1 = 0.064 iC3H7OH(IP) = (CH3)CO(AC)+H2 K2 = 0.076 iC3H7OH(IP) = CHHO(PR) +H2 K3 =

12、 0.00012 后續(xù)加工步驟需要正丙醇，雖然可含丙酮，但丙醛含量不能超過0.05(mol)%。在上述反應(yīng)條件下，是否存在違反這種規(guī)定的可能性？ 數(shù)學(xué)模型：各反應(yīng)的化學(xué)平衡方程如下,function Cha2demo7 x0 = 0.05 0.2 0.01; x = fsolve(EquiC3,x0); CAC=x(3)/sum(x) if CAC0.05 disp(The AC concentration could not be over 0.05%) else disp(The AC concentration could be over 0.05%) end function f =

13、EquiC3(x) f1 = x(1)-0.064*(1-x(1)-x(2)-x(3); f2 = x(2)*(x(2)+x(3)-0.076*(1-x(1)-x(2)-x(3)*(1+x(2)+x(3); f3 = x(3)*(x(2)+x(3)-0.00012*(1-x(1)-x(2)-x(3)*(1+x(2)+x(3); f = f1 f2 f3; 結(jié)果：The AC concentration could not be over 0.05%,function Cha2demo5 rho = 961;G = 6.67;L = 10;eps = 5e-6;deltaP = 15e3;k =

14、 1.8;n = 0.64;K = 1.48;n1 = n; K1 = K*(3*n+1)/(4*n)n; D1 = (32*K1*L*8(n1-1)*(4*G/pi/rho)n1)(1/(3*n1+1); D2 = 2*D1; delta = 1e-4; while abs(D2 - D1)/D2) delta u = (G/rho)/(pi*D12/4); Regen = D1n1*u(2-n1)*rho/(8(n1-1)*K1); if Regen 2100 f = FrictFactor(Regen,n1); end Le = k*D1/(4*f); D1 = D2; D2 = (2*f*(L+Le)/(rho*deltaP)*(4*G/pi)2)0.2; end fprintf(t管道直徑為D = %.4f %sn,D2,m) fprintf(t摩擦因子為f = %.4f %s,f,m) % - function f = FrictFactor(Regen,n1) f0 = 16/Regen; %4.53 f = fzero(fFunc,f0,Regen,n1) % - function F = fFunc(f) F = 1/sqrt(f) - 4.0