计算机化工应用超越方程求解ppt课件.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 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 NextsElseEndIfNextt 200w 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 NextkNextjNexti100 x 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 5 具体应用 运行光盘中的程序 按提示输入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 6 方法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 7 规划求解器求解线性方程组 8 方法3 matlab求解 matlab的命令窗口输入以下命令 a 1234 3421 2152 6121 b 10101010 x a b或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 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 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日 欢迎读者调用 如有问题请告知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 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