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

.,1,计算机化工应用第二讲,华南理工大学化学与化工学院方利国,.,2,问题1线性方程组求解策略,.,3,方法一,利用主元最大高斯消去法，只要有解，就一定可以计算出来，参考光盘第3章程序，输入对应参数即可。,.,4,PrivateSubCommand1_Click()Dimi,j,m,nAsIntegerDima(),z(),x(),w,aa(),s,t,k,ln=InputBox(n)ReDima(n+2,2*n),z(n+2,2*n),x(n+1),aa(n+2,2*n)Fori=1TonForj=1Ton+1a(i,j)=InputBox(输入系数矩阵A(3421;2152;6121;b=10101010;x=ab或x=inv(a)*b计算机输出以下结果：x=1.00001.00001.00001.0000,.,9,问题2非线性方程求解,对于非线性方程或方程组的求解方法和线性方程组求解一样，也有各种迭代方法，如直接迭代、松弛迭代，同时由于非线性方程还可以利用方程本身的信息，延伸出各种加速迭代的方法，如牛顿（newton）迭代、布罗伊登（broyden）迭代、割线（secant）迭代、韦格斯坦（wegsten）迭代。这些迭代法既是非线性方程（组）的求解方法，其实也是大型化工流程模拟时的收敛策略，读者在使用aspenplus等软件模拟时，如果使用软件默认的收敛策略无法收敛时，可以考虑改变收敛方法。,.,10,方法1二分法,.,11,实际问题求解,.,12,Dima2,a3,a4,a5,b2,b3,b5,bvAsDoubleDimc2,c3,c5,tc,t,pAsDoublePrivateSubCommand1_Click()Dima,b,x,x1,x2,y,k,y1,y2AsDoublea2=-4391473.1a3=233734790a4=-8196792900a5=113229830000b2=4501.7239b3=-102972.05b5=74758927bv=20.101853c2=-60767617c3=5081973600c5=-3229376000000tc=304.2温度可修改t=InputBox(TEMPRESURE,)t=273.15+t压力可修改p=InputBox(PRESURE)a=bv+1%二分法初始左边起点，保证解的有效性，不同方程，有不同设置b=a,Dob=b+10%二分法右边点，增加的10可以改变。y1=f(a)%利用自定定义函数调用，将x=a代入自定义函数计算y2=f(b)LoopUntily1*y20%找到有解的范围a,b开始计算Dox=(a+b)/2y=f(x)Ify*y10Thenb=xy2=yElsea=xy1=yEndIfLoopUntilAbs(y)0.0000001PrintV(;p;,;t;)=;xEndSubPublicFunctionf(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)EndFunction,.,13,计算结果（程序见第二章）,利用程序计算温度150，压力为1atm时，CO2气体的摩尔体积为34670.52ml/mol，文献中的数值是34669ml/mol，相对误差为0.00438%；计算温度150，压力为300atm时，CO2气体的摩尔体积为86.54ml/mol，文献中的数值是86.334ml/mol，相对误差为0.239%,.,14,方法2-利用excel求解,.,15,方程设置,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,.,16,宏编程,SubMacro1()Fori=1To12Cells(4,5)=Cells(4+i,5)Cells(4,6)=Cells(4+i,6)Range(I7).SelectRange(I4).GoalSeekGoal:=0,ChangingCell:=Range(G4)Range(J10).SelectCells(4+i,7)=Cells(4,7)Cells(4+i,9)=Cells(4,9)NextEndSub,.,17,方法3-matlab求解,functionqiugen%三种方法求单变量方程根,在7.0版本上调试通过%由华南理工大学方利国编写，2011年11月9日%欢迎读者调用，如有问题请告知lgfangclearallclcx0=-4;x1=fsolve(f,x0)x0=1.2;x2=fzero(f,x0)x0=2.2;x3=fzero(f,x0)%下面的系数向量需和自定义方程中的系数一致，注意系数按降幂排列c=11-2515;x3=roots(c)%-functionf=f(x)%读者可以改变下面的方程次数及系数，只要按范例中的模式书写即可f=x3+x2-25*x+15;,.,18,Matlab求解摩尔体积核心代码,forp=1:30(具体程序见分析与合成第3代码）i=i+1;v(i)=fsolve(f,x0,t,p)eer(i)=f(v(i),t,p)endyp=1:1:30plot(yp,v,r,yp,eer,g),xlabel(P(atm),ylabel(V(ml/mol)holdon；gridonfunctionf=f(x,t,p)globala2a3a4a5b2b3b4b5bvc2c3c5tcf=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),.,19,计算结果,.,20,问题3-非线性方程组求解,一般不建议自己编程，可以使用的方法有：1、excel规则求解，其中有多种表达方式2、matlab求解：fsolve()进行求解，也可以用fminsearch()。fminsearch()求解的策略是将非线性方程组改写成fi(X)=0，在构建目标函数J=(fi2)0.5，通过求目标函数最小值得方法来得到非线性方程组的解。,.,21,实际问题-反应平衡计算,已知进料甲烷为1mol,水蒸汽为5mol,反应后总压P=1atm,反应平衡常数为:,.,22,实际问题-反应平衡计算,解:设进料甲烷为a摩尔，反应平衡时有x摩尔甲烷转化成CO,同时生成的CO中又有y摩尔转化成CO2,假设为理想气体，则反应平衡有以下方程成立：,.,23,Matlab程序,a=1%甲烷进料量(具体程序见分析与合成第3代码）p0=0.80.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)fori=1:56%进行56轮的方程求解a=0.5+(i-1)*0.1p0=0.50.2;%尽量保证初值较理想p1=fsolve(f,p0,a);x=p1(1);y=p1(2);F(i)=(3*x+y)/(5+a+2*x);%氢气摩尔分率end,.,24,Matlab程序,yp=0.5:0.1:6plot(yp,F,r),xlabel(甲烷(mol),ylabel(氢气摩尔分率（%)holdongridonFmax=max(F(:);%求数列中的最大值index=find(F(:)=Fmax);%确定最大值所处的位置max_a=0.5+(index-1)*0.1;fprintf(氢气摩尔分率最大时的甲烷量:),disp(max_a)fprintf(氢气摩尔分率最大值:),disp(Fmax)%-functionf=f(p,a)x=p(1);y=p(2);f(1)=(x-y)*(3*