计算机化工应用超越方程求解PPT课件.ppt

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,.,*,计算机化工应用,第二讲,华南理工大学化学与化工学院方利国,1,.,问题,1,线性方程组求解策略,2,.,方法一,利用主元最大高斯消去法，只要有解，就一定可以计算出来，参考光盘第,3,章程序，输入对应参数即可。,3,.,Private Sub Command1_Click(),Dim i,j,m,n As Integer,Dim a(),z(),x(),w,aa(),s,t,k,l,n=InputBox(n),ReDim a(n+2,2*n),z(n+2,2*n),x(n+1),aa(n+2,2*n),For i=1 To n,For j=1 To n+1,a(i,j)=InputBox(,输入系数矩阵,A(&i&,&j&),Next j,Next I,For i=1 To n,If i=n Then GoTo 200,For t=i+1 To n,If Abs(a(i,i)a=1 2 3 4;3 4 2 1;2 1 5 2;6 1 2 1;,b=10 10 10 10;,x=ab,或,x=inv(a)*b,计算机输出以下结果：,x=,1.0000,1.0000,1.0000,1.0000,8,.,问题,2,非线性方程求解,对于非线性方程或方程组的求解方法和线性方程组求解一样，也有各种迭代方法，如直接迭代、松弛迭代，同时由于非线性方程还可以利用方程本身的信息，延伸出各种加速迭代的方法，如牛顿（,newton,）迭代、布罗伊登（,broyden,）迭代、割线（,secant,）迭代、韦格斯坦（,wegsten,）迭代。这些迭代法既是非线性方程（组）的求解方法，其实也是大型化工流程模拟时的收敛策略，读者在使用,aspen plus,等软件模拟时，如果使用软件默认的收敛策略无法收敛时，可以考虑改变收敛方法。,9,.,方法,1,二分法,10,.,实际问题求解,11,.,Dim a2,a3,a4,a5,b2,b3,b5,bv As Double,Dim c2,c3,c5,tc,t,p As Double,Private Sub Command1_Click(),Dim a,b,x,x1,x2,y,k,y1,y2 As Double,a2=-4391473.1,a3=233734790,a4=-8196792900,a5=113229830000,b2=4501.7239,b3=-102972.05,b5=74758927,bv=20.101853,c2=-60767617,c3=5081973600,c5=-3229376000000,tc=304.2,温度可修改,t=InputBox(TEMPRESURE,),t=273.15+t,压力可修改,p=InputBox(PRESURE),a=bv+1%,二分法初始左边起点，保证解的有效性，不同方程，有不同设置,b=a,Do,b=b+10%,二分法右边点，增加的,10,可以改变。,y1=f(a)%,利用自定定义函数调用，将,x=a,代入自定义函数计算,y2=f(b),Loop Until y1*y2 0%,找到有解的范围,a,b,开始计算,Do,x=(a+b)/2,y=f(x),If y*y1 0 Then,b=x,y2=y,Else,a=x,y1=y,End If,Loop Until Abs(y)0.0000001,Print V(;p;,;t;)=;x,End Sub,Public Function f(x)%,自定义函数，读者需根据不同情况改写,f=p-(82.06*t/(x-bv)+(a2+b2*t+c2*Exp(-5*t/tc)/(x-bv)2+(a3+b3*t+c3*Exp(-5*t/tc)/(x-bv)3),f=f-(a4/(x-bv)4+(a5+b5*t+c5*Exp(-5*t/tc)/(x-bv)5),End Function,12,.,计算结果（程序见第二章）,利用程序计算温度,150,，压力为,1atm,时，,CO,2,气体的摩尔体积为,34670.52ml/mol,，文献中的数值是,34669 ml/mol,，相对误差为,0.00438%,；,计算温度,150,，压力为,300atm,时，,CO,2,气体的摩尔体积为,86.54ml/mol,，文献中的数值是,86.334 ml/mol,，相对误差为,0.239%,13,.,方法,2-,利用,excel,求解,14,.,方程设置,I4=F4-C16*E4/(G4-C11),-(C4+C8*E4+C12*EXP(-5*E4/C15)/(G4-C11)2,-(C5+C9*E4+C13*EXP(-5*E4/C15)/(G4-C11)3,-C6/(G4-C11)4,-(C7+C10*E4+C14*EXP(-5*E4/C15)/(G4-C11)5,15,.,宏编程,Sub Macro1(),For i=1 To 12,Cells(4,5)=Cells(4+i,5),Cells(4,6)=Cells(4+i,6),Range(I7).Select,Range(I4).GoalSeek Goal:=0,ChangingCell:=Range(G4),Range(J10).Select,Cells(4+i,7)=Cells(4,7),Cells(4+i,9)=Cells(4,9),Next,End Sub,16,.,方法,3-matlab,求解,function qiugen,%,三种方法求单变量方程根,在,7.0,版本上调试通过,%,由华南理工大学方利国编写，,2011,年,11,月,9,日,%,欢迎读者调用，如有问题请告知,lgfang,clear all,clc,x0=-4;,x1=fsolve(f,x0),x0=1.2;,x2=fzero(f,x0),x0=2.2;,x3=fzero(f,x0),%,下面的系数向量需和自定义方程中的系数一致，注意系数按降幂排列,c=1 1-25 15;,x3=roots(c),%-,function f=f(x),%,读者可以改变下面的方程次数及系数，只要按范例中的模式书写即可,f=x3+x2-25*x+15;,17,.,Matlab,求解摩尔体积核心代码,for p=1:30(,具体程序见分析与合成第,3,代码）,i=i+1;,v(i)=fsolve(f,x0,t,p),eer(i)=f(v(i),t,p),end,yp=1:1:30,plot(yp,v,r,yp,eer,g),xlabel(P(atm),ylabel(V(ml/mol),hold on,；,grid on,function f=f(x,t,p),global a2 a3 a4 a5 b2 b3 b4 b5 bv c2 c3 c5 tc,f=p-(82.06*t/(x-bv)-(a2+b2*t+c2*exp(-5*t/tc)/(x-bv)2+(a3+b3*t+c3*exp(-5*t/tc)/(x-bv)3)-(a4/(x-bv)4+(a5+b5*t+c5*exp(-5*t/tc)/(x-bv)5),18,.,计算结果,19,.,问题,3-,非线性方程组求解,一般不建议自己编程，可以使用的方法有：,1,、,excel,规则求解，其中有多种表达方式,2,、,matlab,求解：,fsolve(),进行求解，也可以用,fminsearch(),。,fminsearch(),求解的策略是将非线性方程组改写成,f,i,(,X,)=0,，在构建目标函数,J=(,f,i,2,),0.5,，通过求目标函数最小值得方法来得到非线性方程组的解。,20,.,实际问题,-,反应平衡计算,已知进料甲烷为,1mol,水蒸汽为,5mol,反应后总压,P=1atm,反应平衡常数为,:,21,.,实际问题,-,反应平衡计算,解,:,设进料甲烷为,a,摩尔，反应平衡时有,x,摩尔甲烷转化成,CO,同时生成的,CO,中又有,y,摩尔转化成,CO,2,假设为理想气体，则反应平衡有以下方程成立：,22,.,Matlab,程序,a=1%,甲烷进料量,(,具体程序见分析与合成第,3,代码,）,p0=0.8 0.2%,给定,x,y,的初值，初值很重要，如果偏差太大，可能不收敛，需根据物理意义给定,p=fsolve(f,p0,a);%a,作为参数传递进入自定义函数,x=p(1);,y=p(2);,MF(1)=(a-x)/(5+a+2*x);%,甲烷摩尔分率,MF(2)=(5-x-y)/(5+a+2*x);%,水摩尔分率,MF(3)=(x-y)/(5+a+2*x);%,一氧化碳摩尔分率,MF(4)=y/(5+a+2*x);%,二氧化碳摩尔分率,MF(5)=(3*x+y)/(5+a+2*x);%,氢气摩尔分率,fprintf(x,及,y,值,:),disp(p),fprintf(,摩尔分率,:),disp(MF),for i=1:56%,进行,56,轮的方程求解,a=0.5+(i-1)*0.1,p0=0.5 0.2;%,尽量保证初值较理想,p1=fsolve(f,p0,a);,x=p1(1);,y=p1(2);,F(i)=(3*x+y)/(5+a+2*x);%,氢气摩尔分率,end,23,.,Matlab,程序,yp=0.5:0.1:6,plot(yp,F,r),xlabel(,甲烷,(mol),ylabel(,氢气摩尔分率（,%),hold on,grid on,Fmax=max(F(:);%,求数列中的最大值,index=find(F(:)=Fmax);%,确定最大值所处的位置,max_a=0.5+(index-1)*0.1;,fprintf(,氢气摩尔分率最大时的甲烷量,:),disp(max_a),fprintf(,氢气摩尔分率最大值,:),disp(Fmax),%-,function f=f(p,a),x=p(1);,y=p(2);,f(1)=(x-y)*(3*x+y)3/(a-x)*(5-x-y)*(5+a+2*x)2)-0.9618;,f(2)=y*(3*x+y)/(x-y)*(5-x-y)-2.7;,f=f(1)f(2);,24,.,计算结果,25,.,Excel,计算,26,.,Excel,计算,a,x,y,f1,f2,j,3.3,2.030271,0.901484,-1.8E-06,5.58E-07,3.68E-12,0,0,CH4,H2O,CO,CO2,H2,10.27244,16.73264,9.132183,7.293239,56.56951,100,=(D4-E4)*(3*D4+E4)3/(C4-D4)*(5-D4-E4)*(5+C4+2*D4)2)-0.9618,27,.,Excel,计算,28,.,