计算机化工应用超越方程求解备课讲稿

计算机化工应用超越方程求解方法一利用主元最大高斯消去法，只要有解，就一定可以计算出来，参考光盘第3章程序，输入对应参数即可。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("&i&","&j&")")NextjNextIFori=1TonIfi=nThenGoTo200Fort=i+1TonIfAbs(a(i,i))<Abs(a(t,i))ThenFors=iTon+1aa(t,s)=a(i,s)a(i,s)=a(t,s)a(t,s)=aa(t,s)NextsElseEndIfNextt200w=a(i,i)Forj=1Ton+1a(i,j)=a(i,j)/wNextjIfi=nThenGoTo100Forj=i+1TonFork=i+1Ton+1z(i,k)=a(i,k)*a(j,i)a(j,k)=a(j,k)-z(i,k)NextkNextjNexti100x(n+1)=0Fork=nTo1Step-1s=0Forj=k+1Tons=s+a(k,j)*x(j)Nextjx(k)=a(k,n+1)-sPrint"x(";k;")=";x(k)NextkEndSub具体应用运行光盘中的程序，按提示输入4；1，2，3，4，10；3，4，2，1，10；2，1，5，2，10；6，1，2，1，10。系统计算得到：

x1=1，x2=1，x3=1，x4=1。

方法2-excel规划求解

L4=$B$14*B4+$C$14*C4+$D$14*D4+$E$14*E4+$F$14*F4+$G$14*G4+$H$14*H4+$I$14*I4+$J$14*J4+$K$14*K4规划求解器求解线性方程组方法3-matlab求解matlab的命令窗口输入以下命令：>>a=[1234;3421;2152;6121];>>b=[10101010]';>>x=a\b或x=inv(a)*b计算机输出以下结果：x=1.00001.00001.00001.0000问题2非线性方程求解

对于非线性方程或方程组的求解方法和线性方程组求解一样，也有各种迭代方法，如直接迭代、松弛迭代，同时由于非线性方程还可以利用方程本身的信息，延伸出各种加速迭代的方法，如牛顿（newton）迭代、布罗伊登（broyden）迭代、割线（secant）迭代、韦格斯坦（wegsten）迭代。这些迭代法既是非线性方程（组）的求解方法，其实也是大型化工流程模拟时的收敛策略，读者在使用aspenplus等软件模拟时，如果使用软件默认的收敛策略无法收敛时，可以考虑改变收敛方法。方法1—二分法实际问题求解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*y2<0%找到有解的范围[a,b]'开始计算

Dox=(a+b)/2y=f(x)Ify*y1<0Thenb=xy2=yElsea=xy1=yEndIfLoopUntilAbs(y)<0.0000001Print"V(";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计算结果（程序见第二章）利用程序计算温度150℃，压力为1atm时，CO2气体的摩尔体积为34670.52ml/mol，文献中的数值是34669ml/mol，相对误差为0.00438%；计算温度150℃，压力为300atm时，CO2气体的摩尔体积为86.54ml/mol，文献中的数值是86.334ml/mol，相对误差为0.239%方法2-利用excel求解方程设置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宏编程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方法3-matlab求解functionqiugen%三种方法求单变量方程根,在7.0版本上调试通过%由华南理工大学方利国编写，2011年11月9日%欢迎读者调用，如有问题请告知lgfang@clearallclcx0=-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=x^3+x^2-25*x+15;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)计算结果问题3-非线性方程组求解

一般不建议自己编程，可以使用的方法有：1、excel规则求解，其中有多种表达方式2、matlab求解：fsolve()进行求解，也可以用fminsearch()。fminsearch()求解的策略是将非线性方程组改写成fi(X)=0，在构建目标函数J=(∑fi2)0.5，通过求目标函数最小值得方法来得到非线性方程组的解。实际问题-反应平衡计算已知进料甲烷为1mol,水蒸汽为5mol,反应后总压P=1atm,反应平衡常数为:实际问题-反应平衡计算解:设进料甲烷为a摩尔，反应平衡时有x摩尔甲烷转化成CO,同时生成的CO中又有y摩尔转化成CO2,假设为理想气体，则反应平衡有以下方程成立：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);%氢气摩尔分率endMatlab程序yp=0.5:0.1:6plot(yp,F,'r'),xlabel('甲烷(mol)'),ylabel('氢气摩尔分率（%)')holdongridonFmax=max(F(:));%求数列中