Matlab与化学化工计算.pptVIP

计算机在化学化工中的应用七Matlab与化学化工计算 本节要点 本章背景Matlab基础方程组求解数据插值作业 问题的提出 MATLAB语言与其它语言的关系仿佛和C语言与汇编语言的关系一样计算机语言的发展标志着计算机语言向 智能化 方向发展 被称为第四代编程语言 1Matlab基础知识 1 1Matlab简介 1967年由ClereMaler用FORTRAN语言设计和编写1984年Mathworks公司用C语言完成了Matlab的商业化版本并推向市场经过20余年的改进 Matlab已发展成为一个具有极高通用性的 带有众多实用工具的运算平台 成为国际上广泛认可的优秀科学计算软件 Matlab的发展 1984年 MATLAB第1版 DOS版 1992年 MATLAB4 0版1994年 MATLAB4 2版1997年 MATLAB5 0版1999年 MATLAB5 3版2000年 MATLAB6 0版2001年 MATLAB6 1版2002年 MATLAB6 5版2004年 MATLAB7 0版 告别DOS版 1993年MathWorks公司从加拿大滑铁卢大学购得Maple的使用权 推出了符号计算工具包 5 0的MATLAB拥有更丰富的数据类型和结构 更友善的面向对象 更加快速精良的图形可视 更广博的数学和数据分析资源 更多的应用开发工具 Matlab的优点 语法简单易学 编程效率高高质量 高可靠的数值计算能力强大的矩阵运算能力高级图形和数据可视化处理能力提供600多个常用算法内建函数 以及众多面向应用的工具箱 Matlab二维作图 Matlab三维作图 1 2Matlab的界面 1 3Matlab的帮助功能 联机帮助系统命令窗口查询helplookfor联机演示系统Demos Help 下拉菜单中 FullProductFamilyHelp 命令打开联机帮助系统 若不知函数确切名 可 Lookfor关键词 可查 help Help全部主题 Help指定函数 例7 1 查找包含 diff 关键词的函数 lookfordiffSETDIFFSetdifference DIFFDifferenceandapproximatederivative POLYDERDifferentiatepolynomial DDE23Solvedelaydifferentialequations DDEs withconstantdelays DDESDSolvedelaydifferentialequations DDEs withgeneraldelays DEVALEvaluatethesolutionofadifferentialequationproblem 用户输入的命令 查询结果 2线性方程组求解 2 1线性方程组的一般形式 在应用中 常常把线性方程组写成AX b的一般形式 其中 2 2线性方程组解的判断 齐次线性方程组AX 0 其解的情况可以通过系数矩阵A的秩和未知数个数n的关系来判断如果系数矩阵的秩为n 方程组只有零解 x 0如果系数矩阵的秩小于n 方程组有无穷多解如果系数矩阵的秩大于n 方程组无解 非其次线性方程组解的情况 非齐次线性方程组AX b 根据系数矩阵A的秩 增广矩阵B Ab 的秩和未知数个数n的关系来判断其解的情况如果系数矩阵A的秩等于增广矩阵B的秩且等于n 方程组有唯一解如果系数矩阵A的秩等于增广矩阵B的秩且小于n 方程组有无穷多解如果系数矩阵A的秩小于增广矩阵B的秩 方程组无解 例7 2判断方程解的情况 解 在Matlab中输入 a 1 24 211 11 1 rank a ans 2齐次线性方程组系数矩阵A的秩为2 小于未知数个数3 方程组有无穷多解 计算系数矩阵A的秩 不能少 例7 2 2 解 a 7028 0281 280196 b 1 39 7 b为列向量 故输入行向量后转置 rank a 计算系数矩阵A的秩ans 3 rank ab 计算增广矩阵 Ab 的秩ans 3非齐次线性方程组系数矩阵A的秩为3 增广矩阵的秩为3 等于未知数个数3 方程组有唯一解 是Matlab的注释符 后的语句作为注释处理 2 3线性方程组直接求解 例7 3求以下方程组的解步骤b1矩阵除法 验证解a判断解的情况b2逆矩阵法b3rref 例7 4 求下列方程组的解 视频演示 3数据插值 3 1数据插值简介 在工程领域 许多实验数据常以列表函数或表格的形式存在 如水黏度随温度的列表函数在实际使用时 有时需要获得介于表中两个温度结点之间 如15 25 的黏度值 而这些数据未在表中出现 需要我们根据已知的数据估算出表中未出现的温度点的黏度数值 这一技术称为插值技术 插值的数学定义 已知由g X 可能未知或非常复杂 产生的n 1个离散数据 xi yi i 0 1 2 n 且这n 1个互异插值结点满足a x0 x1 x2 xn b 在插值区间 a b 内寻找一个相对简单的函数f x 使其满足插值条件f xi yi i 0 1 2 n 再利用已求得的f x 计算任一非插值结点x 处的近似值y f x 其中f x 称为插值函数 g x 称为被插值函数从计算的观点看 插值就是用一个简单函数在某种误差范围内近似的代替原目标函数关系式 3 2插值方法 线性插值二次插值其他插值方法最近 nearest 插值法样条曲线 spline 法埃尔米特 Hermite 法 3 2 1线性插值 又称两点插值已知两个数据点 x0 y0 x1 y1 x0 x1 求对应于x x0 x x1 的y值解法 由x0 y0 x1 y1构造直线方程并求取在该点的函数值 线性插值的优点是简单 快捷 特别是对于插值结点间距较小的情况可以取得令人满意的精度 3 2 2二次插值 又称拉格朗日三点差值根据三个已知点 x0 y0 x1 y1 x2 y2 x0 x1 x2 构造二次多项式插值函数y a0 a1x a2x2 并用该函数计算在x处的y值二次插值公式 3 3 1使用Matlab进行数据插值 一维插值只有一个自变量的插值Matlab提供的一维插值函数是interp1常用语法 YI interp1 X Y XI method 式中X Y 为已知数据点的x y值 XI 为待插值数据点的x值 YI 为返回的插值结果 method 用于指定所采用的插值方法 插值函数interp1提供的插值方法 不同方法插值的结果 在 0 2 区间内生成11个等距的离散点 计算函数y sin x 的数值 分段三次Hermite插值 分段三次样条插值 线性插值 最近插值 例7 5 用函数y ex生成以下离散数据 使用Matlab的不同插值方法计算x 2 552 632 772 86 处的函数值 并与真实值进行比较 视频演示 例7 5计算结果的比较 最接近真实值 例7 6 已知水在20 21 22 23 的饱和蒸汽压分别为17 54 18 65 19 83 21 07mmHg 求20 5 21 5 22 5和24 时水的饱和蒸汽压各是多少 已知24 时水的饱和蒸汽压为22 38mmHg 视频演示 3 3 2多维插值 具有多个自变量的插值二维插值三维插值高维插值 插值函数 Method选项 例7 7 函数z ex sin y y 1生成表7 6中的离散数据 应用不同插值方法计算在x 0 36 y 1 9处的z值 并与真值作比较 例7 7插值结果 视频演示 最接近真实值 4非线性方程的求解 4 1非线性方程数值求解 常见的非线性方程 组 有两种形式或 直接迭代法 韦斯顿迭代法求解 牛顿迭代法求解 4 1 1直接迭代法 先设定x的初值x x0 代入方程计算x1 f x0 再把x1作为新的初值代入方程计算x2 f x1 直至求得的xn 1与xn足够接近 称为收敛 xn即为方程的根直接迭代法求解可写为如下形式 直接迭代法的优点是形式简单 易于编程实现 缺点是计算量大 收敛速度慢 一般可通过改进初值 降低收敛要求等方法提高其收敛速度 也可采用其它方法进行求解 收敛判断准则 绝对偏差相对偏差半相对偏差 是用户指定的一个很小的正数 确定适当的 取值有一定难度 优点在于 的选取不受方程根的数值大小的影响 一般取 0 001 是判断收敛的一个较好的方法 的取值一般为0 0001 0 001 4 1 2韦格斯顿迭代法 韦斯顿法对直接迭代做了改进 使用前两个计算点的信息进行求解 其迭代形式为韦格斯顿法迭代时需要前面两个计算点的数据 可先执行一次直接迭代法计算获得 4 1 3牛顿迭代法 又称切线法若将f x 在其根附近进行泰勒级数展开 并取级数的线性部分作为f x 的近似值 可得由上式可得牛顿迭代公式如函数的一阶导数难以求得 可用差商作为近似导数 4 2用Matlab求解非线性方程 Matlab求解非线性方程 组 的函数fzero 一元非线性方程的求解fsolve 用于非线性方程组的求解fzero函数用法x fzero fun x0 fsolve函数用法x fsolve fun x0 fun 单变量实值函数 可以是Matlab内部函数或用户自定义函数x0 若x0是一个单个的数值 系统会将其作为求解的初值 在其附近寻找解 若x0是一个二维向量 且fun x0 1 和fun x0 2 符号相反 Matlab将会在x0 1 和x0 2 区间内寻找零点 fun 用户自定义函数 返回给定变量x时方程 组 的值y fun x x0 初值矩阵 对于fzero和fsolve函数 给定适当的初值对问题的求解至关重要 若初值选择不当 将无法得到正确的解 一般可根据经验或简化计算获得合适的初值 例7 8 试用维里方程计算200 1 013MPa的异丙醇蒸汽的摩尔体积V与压缩因子Z 已知异丙醇的维里系数实验值B 388cm3 mol 1 C 26000cm6 mol 2 视频演示 例7 9 600K下由CH3Cl和H2O反应生成CH3OH 存在下列平衡CH3Cl g H2O g CH3OH g HCl g 1 2CH3OH g CH3 2O g H2O g 2 已知该温度下Kp 1 0 00154 Kp 2 10 6 今以等摩尔的CH3Cl g 和H2O g 开始反应 求CH3Cl的平衡转化率 视频演示 5常微分方程 组 求解 5 1化工中的常微分方程 组 微分方程中只有一个自变量的方程称为常微分方程 自变量个数为两个或两个以上的微分方程称为偏微分方程常微分方程初值问题 给定微分方程及初值条件边值问题 给定微分方程及边界条件常微分方程 组 的解法解析法和数值法 常用 初值问题 记为或 5 2常微分方程 组 数值解法 欧拉公式梯形公式龙格 库塔法常微分方程组的数值解法 5 2 1欧拉公式 若常微分方程初值问题的求解区间为 将其等分为m步 步长 记 相应xn处的函数值为yn 则yn可由下式计算向前欧拉公式向后欧拉公式中心欧拉公式 若yn 1同时出现在等号的两侧 称为隐式欧拉公式 无法直接求解 一般需采用迭代法计算 5 2 2梯形公式 梯形公式也是隐式格式 需要迭代求解先用显式公式算出初值 再用隐式公式进行一次或数次修正 这一过程称为预估 校正过程公式为合并为 5 2 3龙格 库塔法 工程应用中求解常微分方程最常用的一种有效方法 其计算精度和运算速度较快 易于编程 常用的有二阶 三阶 四阶龙格 库塔公式二阶龙格 库塔公式 常见形式或 三阶龙格 库塔公式 常见形式 四阶龙格 库塔公式 常见形式 5 2 4常微分方程组的数值解法 将由m个一阶方程组成的常微分初值问题其中 可由前边所述的解常微分方程的各个方法求解 写为向量形式 5 2 5高阶常微分方程数值解法 可把高阶常微分方程转化为一阶常微分方程组求解 例如三阶常微分方程令将三阶方程化为一阶方程 5 3用Matlab求解常微分方程 例7 10 在间歇反应器中进行液相反应制备产物B 其反应网络如图7 7所示 反应温度为224 6 反应物X大量过剩 各反应均为一级动力学关系 各步反应的k0i Eai见表 试给出0 10000秒各产物的浓度变化规律 初始条件为 t 0 CA 1kmol m3 CB CC CD CE 0kmol m3 反应网络图 例7 10条件图 反应网络图 参数取值