计算机化工应用超越方程求解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) Abs(a(t, i) Then For s = i To n + 1 aa(t, s) = a(i, s) a(i, s) = a(t, s) a(t, s) = aa(t, s) Next s Else End If Next t,200 w = a(i, i) For j = 1 To n + 1 a(i, j) = a(i, j) / w Next j If i = n Then GoTo 100 For j = i + 1 To n For k = i + 1 To n + 1 z(i, k) = a(i, k) * a(j, i) a(j, k) = a(j, k) - z(i, k) Next k Next j Next i 100 x(n + 1) = 0 For k = n To 1 Step -1 s = 0 For j = k + 1 To n s = s + a(k, j) * x(j) Next j x(k) = a(k, n + 1) - s Print “x(“; k; “)=“; x(k) Next k End Sub,4,具体应用,运行光盘中的程序，按提示输入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。,5,方法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,6,规划求解器求解线性方程组,7,方法3-matlab求解,mat lab的命令窗口输入以下命令： 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时，CO2气体的摩尔体积为34670.52ml/mol，文献中的数值是34669 ml/mol，相对误差为0.00438%；计算温度150，压力为300atm时，CO2气体的摩尔体积为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()求解的策略是将非线性方程组改写成fi(X)=0，在构建目标函数J=(fi2)0.5，通过求目标函数最小值得方法来得到非线性方程组的解。,20,实际问题-反应平衡计算,已知进料甲烷为1mol,水蒸汽为5mol,反应后总压P=1atm,反应平衡常数为:,21,实际问题-反应平衡计算,解:设进料甲烷为a 摩尔，反应平衡时有x摩尔甲烷转化成CO,同时生成的CO中又有y摩尔转化成CO2, 假设为理想气体，则反应平衡有以下方程成立：,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计算,=(D4-E4)*(3*D4+E4)3/(C4-D4)*(5-D4-E4)*(5+C4+2*D4)2)-0.9618,27,Excel计算,28,