Matlab与化学化工计算

1、,计算机在化学化工中的应用 七 matlab与化学化工计算,本节要点,本章背景 matlab基础 方程组求解 数据插值 作业,1 matlab 基础知识,1.1 matlab 简介,1967年由clere maler用fortran语言设计和编写 1984年mathworks公司用c语言完成了matlab的商业化版本并推向市场 经过20余年的改进，matlab已发展成为一个具有极高通用性的、带有众多实用工具的运算平台，成为国际上广泛认可的优秀科学计算软件,matlab 的发展,1984年，matlab第1版(dos版) 1992年，matlab 4.0版 1994年，matlab 4.2版 1

2、997年，matlab 5.0版 1999年，matlab 5.3版 2000年，matlab 6.0版 2001年，matlab 6.1版 2002年，matlab 6.5版 2004年，matlab 7.0版,告别dos版,1993年mathworks公司从加拿大滑铁卢大学购得maple的使用权，推出了符号计算工具包,5.0的matlab拥有更丰富的数据类型和结构、更友善的面向对象、更加快速精良的图形可视、更广博的数学和数据分析资源、更多的应用开发工具,matlab三维作图,1.2 matlab 的界面,1.3 matlab 的帮助功能,联机帮助系统 命令窗口查询 help lookfor

3、 联机演示系统 demos,“help”下拉菜单中“full product family help”命令打开联机帮助系统,若不知函数确切名，可“lookfor关键词”可查,help,help全部主题,help指定函数,例7-1,查找包含“diff”关键词的函数 lookfor diff setdiff set difference. diff difference and approximate derivative. polyder differentiate polynomial. dde23 solve delay differential equations (ddes) with

4、constant delays. ddesd solve delay differential equations (ddes) with general delays. deval evaluate the solution of a differential equation problem. ,用户输入的命令,查询结果,2 线性方程组求解,2.1 线性方程组的一般形式,在应用中，常常把线性方程组 写成ax=b的一般形式，其中,2.2 线性方程组解的判断,齐次线性方程组ax=0，其解的情况可以通过系数矩阵a的秩和未知数个数n的关系来判断 如果系数矩阵的秩为n，方程组只有零解，x=0 如果系

5、数矩阵的秩小于n，方程组有无穷多解 如果系数矩阵的秩大于n，方程组无解,非其次线性方程组解的情况,非齐次线性方程组ax=b，根据系数矩阵a的秩、增广矩阵b=a b的秩和未知数个数n的关系来判断其解的情况 如果系数矩阵a的秩等于增广矩阵b的秩且等于n ，方程组有唯一解 如果系数矩阵a的秩等于增广矩阵b的秩且小于n ，方程组有无穷多解 如果系数矩阵a的秩小于增广矩阵b的秩，方程组无解,例7-2 判断方程解的情况,解：在matlab中输入 a=-1 -2 4; 2 1 1; 1 1 -1; rank(a) ans = 2 齐次线性方程组系数矩阵a的秩为2，小于未知数个数3，方程组有无穷多解,计算系数

6、矩阵a的秩,;不能少,例7-2（2）,解： a=7 0 28; 0 28 1; 28 0 196; b=1 -39 -7; %b为列向量，故输入行向量后转置 rank(a) %计算系数矩阵a的秩 ans = 3 rank(a b) %计算增广矩阵a b的秩 ans = 3 非齐次线性方程组系数矩阵a的秩为3，增广矩阵的秩为3，等于未知数个数3，方程组有唯一解。,“%”是matlab的注释符，%后的语句作为注释处理,2.3 线性方程组直接求解,例7-3 求以下方程组的解 步骤 b1 矩阵除法，验证解 a 判断解的情况 b2 逆矩阵法 b3 rref,例7-4,求下列方程组的解,视频演示,3 数据

7、插值,3.1 数据插值简介,在工程领域，许多实验数据常以列表函数或表格的形式存在，如水黏度随温度的列表函数 在实际使用时，有时需要获得介于表中两个温度结点之间（如15，25）的黏度值。而这些数据未在表中出现，需要我们根据已知的数据估算出表中未出现的温度点的黏度数值，这一技术称为插值技术,插值的数学定义,已知由g(x) (可能未知或非常复杂)产生的n+1个离散数据(xi, yi), i=0,1,2,n，且这n+1个互异插值结点满足a=x0x1x2xn=b，在插值区间a,b内寻找一个相对简单的函数f(x)，使其满足插值条件f(xi)=yi, i=0,1,2,n。再利用已求得的 f(x)计算任一非插

8、值结点x*处的近似值y*=f(x*)。其中f(x)称为插值函数，g(x)称为被插值函数 从计算的观点看，插值就是用一个简单函数在某种误差范围内近似的代替原目标函数关系式,3.2 插值方法,线性插值 二次插值 其他插值方法 最近（nearest）插值法 样条曲线（spline）法 埃尔米特（hermite）法,3.2.1 线性插值,又称两点插值 已知两个数据点x0, y0, x1, y1 (x0x1)，求对应于x (x0xx1)的y值 解法：由x0,y0,x1,y1构造直线方程 并求取在该点的函数值,线性插值的优点是简单，快捷，特别是对于插值结点间距较小的情况可以取得令人满意的精度,3.2.2

9、二次插值,又称拉格朗日三点差值 根据三个已知点 x0, y0, x1,y1 , x2,y2 (x0x1x2)，构造二次多项式插值函数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个

10、等距的离散点，计算函数y=sin(x)的数值,分段三次hermite插值,分段三次样条插值,线性插值,最近插值,例7-5,用函数y=ex生成以下离散数据，使用matlab的不同插值方法计算x=2.55 2.63 2.77 2.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 多维插值,具有多个自变量的插值 二维插值

11、三维插值 高维插值,插值函数,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即为方程的根 直接迭代法求解可写为如下形式,直接迭代法

12、的优点是形式简单，易于编程实现。缺点是计算量大、收敛速度慢。一般可通过改进初值、降低收敛要求等方法提高其收敛速度。也可采用其它方法进行求解,收敛判断准则,绝对偏差 相对偏差 半相对偏差,是用户指定的一个很小的正数，确定适当的取值有一定难度,优点在于的选取不受方程根的数值大小的影响。一般取=0.001,是判断收敛的一个较好的方法，的取值一般为0.00010.001,4.1.2 韦格斯顿迭代法,韦斯顿法对直接迭代做了改进，使用前两个计算点的信息进行求解。其迭代形式为 韦格斯顿法迭代时需要前面两个计算点的数据，可先执行一次直接迭代法计算获得,4.1.3 牛顿迭代法,又称切线法 若将f(x)在其根附近

13、进行泰勒级数展开，并取级数的线性部分作为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)符号相

14、反，matlab将会在x0(1)和x0(2)区间内寻找零点,fun用户自定义函数，返回给定变量x时方程(组)的值y=fun(x) x0 初值矩阵,对于fzero和fsolve函数，给定适当的初值对问题的求解至关重要，若初值选择不当，将无法得到正确的解。一般可根据经验或简化计算获得合适的初值,例7-8,试用维里方程计算200，1.013mpa的异丙醇蒸汽的摩尔体积v与压缩因子z。已知异丙醇的维里系数实验值b=-388cm3mol-1，c=-26000 cm6mol-2,视频演示,例7-9,600k下由ch3cl和h2o反应生成ch3oh，存在下列平衡 ch3cl(g)+h2o(g)=ch3oh(

15、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 常微分方程（组）数值解法,欧拉公式 梯形公式

16、龙格-库塔法 常微分方程组的数值解法,5.2.1 欧拉公式,若常微分方程初值问题的求解区间为 ，将其等分为m步，步长 。记 ，相应xn处的函数值为yn , ，则yn可由下式计算 向前欧拉公式 向后欧拉公式 中心欧拉公式,若yn+1同时出现在等号的两侧，称为隐式欧拉公式，无法直接求解，一般需采用迭代法计算,5.2.2 梯形公式,梯形公式也是隐式格式，需要迭代求解 先用显式公式算出初值，再用隐式公式进行一次或数次修正。这一过程称为预估-校正过程 公式为 合并为,5.2.3 龙格-库塔法,工程应用中求解常微分方程最常用的一种有效方法，其计算精度和运算速度较快，易于编程。常用的有二阶、三阶、四阶龙格-

17、库塔公式 二阶龙格-库塔公式 ，常见形式 或,三阶龙格-库塔公式,常见形式,四阶龙格-库塔公式,常见形式,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=0 kmol/m3,反应网络图,例7-10条件图,反应网络图,参数取值,例7-1