Matlab与化学化工计算.pptx

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,版权所有：河南大学化学化工学院 李谦,2010,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,版权所有：河南大学化学化工学院 李谦,2010,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,版权所有：河南大学化学化工学院 李谦,2010,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,版权所有：河南大学化学化工学院 李谦,2010,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,版权所有：河南大学化学化工学院 李谦,2010,*,计算机在化学化工中应用,七,Matlab,与化学化工计算,1/66,本节关键点,本章背景,Matlab,基础,方程组求解,数据插值,作业,2/66,问题提出,MATLAB,语言与其它语言关系好像和,C,语言与汇编语言关系一样,计算机语言发展,标志着计算机语言向“智能化”方向发展，被称为第四代编程语言,数值运算,解析运算,管理、可视化,智能化,3/66,1 Matlab,基础知识,4/66,1.1 Matlab,介绍,1967,年由,Clere Maler,用,FORTRAN,语言设计和编写,1984,年,Mathworks,企业用,C,语言完成了,Matlab,商业化版本并推向市场,经过20余年改进,，,Matlab已发展成为一个含有极高通用性、带有众多实用工具运算平台，成为国际上广泛认可优异科学计算软件,5/66,Matlab,发展,1984,年，,MATLAB,第,1,版,(DOS,版,)1992,年，,MATLAB 4.0,版,1994,年，,MATLAB 4.2,版,1997,年，,MATLAB 5.0,版1999,年，,MATLAB 5.3,版,年，,MATLAB 6.0,版,年，,MATLAB 6.1,版,年，,MATLAB 6.5,版,年，,MATLAB 7.0,版,告别,DOS,版,1993,年,MathWorks,企业从加拿大滑铁卢大学购得,Maple,使用权，推出了符号计算工具包,5.0,MATLAB,拥有更丰富数据类型和结构、更友善面向对象、愈加紧速精良图形可视、更广博数学和数据分析资源、更多应用开发工具,6/66,Matlab,优点,语法简单易学，编程效率高,高质量、高可靠数值计算能力,强大矩阵运算能力,高级图形和数据可视化处理能力,提供600多个惯用算法内建函数，以及众多面向应用工具箱,7/66,Matlab,二维作图,8/66,Matlab,三维作图,9/66,1.2 Matlab,界面,10/66,1.3 Matlab,帮助功效,联机帮助系统,命令窗口查询,help,lookfor,联机演示系统,Demos,“,Help,”下拉菜单中“,Full Product Family Help,”,命令打开联机帮助系统,若不知函数确切名，可“,Lookfor,关键词”可查,11/66,help,Help,全部主题,Help,指定函数,12/66,例,7-1,查找包含,“diff”,关键词函数,lookfor diff,SETDIFF Set difference.,DIFF Difference and approximate derivative.,POLYDER Differentiate polynomial.,DDE23 Solve delay differential equations(DDEs)with constant delays.,DDESD Solve delay differential equations(DDEs)with general delays.,DEVAL Evaluate the solution of a differential equation problem.,用户输入命令,查询结果,13/66,2,线性方程组求解,14/66,2.1,线性方程组普通形式,在应用中，经常把线性方程组,写成,AX,=,b,普通形式，其中,15/66,2.2,线性方程组解判断,齐次线性方程组,AX,=0,，其解情况能够经过系数矩阵,A,秩和未知数个数,n,关系来判断,假如系数矩阵秩为,n,，方程组只有零解，,x,=0,假如系数矩阵秩小于,n,，方程组有没有穷多解,假如系数矩阵秩大于,n,，方程组无解,16/66,非其次线性方程组解情况,非齐次线性方程组,AX,=,b,，依据系数矩阵,A,秩、增广矩阵,B,=,A,b,秩和未知数个数,n,关系来判断其解情况,假如系数矩阵,A,秩等于增广矩阵,B,秩且等于,n,，方程组有唯一解,假如系数矩阵,A,秩等于增广矩阵,B,秩且小于,n,，方程组有没有穷多解,假如系数矩阵,A,秩小于增广矩阵,B,秩，方程组无解,17/66,例,7-2,判断方程解情况,解：在,Matlab,中输入,a=-1-2 4;2 1 1;1 1-1;,rank(a),ans=,2,齐次线性方程组系数矩阵,A,秩为,2,，小于未知数个数,3,，方程组有没有穷多解,计算系数矩阵,A,秩,;,不能少,18/66,例,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,注释符，,%,后语句作为注释处理,19/66,2.3,线性方程组直接求解,例,7-3,求,以下方程组解,步骤,b,1,矩阵除法，验证解,a,判断解情况,b,2,逆矩阵法,b,3,rref,20/66,例,7-4,求以下方程组解,视频演示,21/66,3,数据插值,22/66,3.1,数据插值介绍,在工程领域，许多试验数据常以列表函数或表格形式存在,，如,水黏度随温度列表函数,在实际使用时，有时需要取得介于表中两个温度结点之间（如,15,，,25,）黏度值。而这些数据未在表中出现，需要我们依据已知数据估算出表中未出现温度点黏度数值，这一技术称为,插值技术,t/,0,10,20,30,40,50,60,/mPas,1.788,1.305,1.004,0.8012,0.6532,0.5492,0.4698,23/66,插值数学定义,已知由,g,(,X,)(,可能未知或非常复杂,),产生,n,+1,个离散数据,(,x,i,y,i,),i=0,1,2,n,，且,这,n,+1,个互异插值结点满足,a,=,x,0,x,1,x,2,x,n,=,b,，在插值区间,a,b,内寻找一个相对简单函数,f,(,x,),，使其满足插值条件,f,(,x,i,)=,y,i,i=0,1,2,n,。再利用已求得,f,(,x,),计算任一非插值结点,x,*,处近似值,y,*,=,f,(,x,*,),。其中,f,(,x,),称为插值函数，,g,(,x,),称为被插值函数,从计算观点看，插值,就,是用一个简单函数在某种误差范围内近似代替原目标函数关系式,24/66,3.2,插值方法,线性插值,二次插值,其它插值方法,最近（,nearest,）插值法,样条曲线（,spline,）法,埃尔米特（,Hermite,）法,25/66,3.2.1,线性插值,又称两点插值,已知两个数据点,x,0,y,0,x,1,y,1,(,x,0,x,1),，求对应于,x,(,x,0,x,x,1,),y,值,解法：由,x,0,y,0,x,1,y,1,结构直线方程,并求取在该点函数值,线性插值优点,是简单，快捷，尤其是对于插值结点间距较小情况能够取得令人满意精度,26/66,3.2.2,二次插值,又称拉格朗日三点差值,依据三个已知点,x,0,y,0,x,1,y,1,x,2,y,2,(,x,0,x,1,x,2,),，结构二次多项式插值函数,y,=,a,0,+,a,1,x,+,a,2,x,2,，并用该函数计算在,x,处,y,值,二次插值公式,27/66,3.3.1,使用,Matlab,进行数据插值,一维插值,只有一个自变量插值,Matlab,提供一维插值函数是,interp1,惯用语法：,YI=interp1,(,X,Y,XI,method,),式中,X,Y,为已知数据点,x,y,值；,XI,为待插值数据点,x,值；,YI,为返回插值结果；,method,用于指定所采取插值方法,28/66,插值函数,interp1,提供插值方法,项目,取值,nearest,最近插值,linear,线性插值，为interp1缺省方法,spline,分段三次样条插值,pchip,分段三次,Hermite,插值,cubic,同,pchip,29/66,不一样方法插值结果,在,0,2,区间内生成,11,个等距离散点，计算函数,y,=sin(,x,),数值,分段三次,Hermite,插值,分段三次样条插值,线性插值,最近插值,30/66,例,7-5,用函数,y,=,e,x,生成以下离散数据，使用,Matlab,不一样插值方法计算,x=2.55 2.63 2.77 2.86,处函数值，并与真实值进行比较,x,2.5,2.6,2.7,2.8,2.9,y,12.1825,13.4637,14.8797,16.4446,18.1741,视频演示,31/66,例,7-5,计算结果比较,插值方法,x,2.55,2.63,2.77,2.86,真实值,12.8071,13.8738,15.9586,17.4615,最近插值,13.4637,13.4637,16.4446,18.1741,线性插值,12.8231,13.8885,15.9752,17.4823,三次样条插值,12.8071,13.8738,15.9586,17.4616,分段三次,Hermite,插值,12.8067,13.8737,15.9588,17.4622,最靠近真实值,32/66,例,7-6,已知水在,20,21,22,23,饱和蒸汽压分别为,17.54,18.65,19.83,21.07mmHg,求,20.5,，,21.5,，,22.5,和,24,时水饱和蒸汽压各是多少？已知,24,时水饱和蒸汽压为,22.38mmHg,视频演示,33/66,3.3.2,多维插值,含有多个自变量插值,二维插值,三维插值,高维插值,函数名,经典应用,说明,interp2,ZI=interp2(X,Y,Z,XI,YI,method),二维插值,interp3,VI=interp3(X,Y,Z,V,XI,YI,ZI,method),三维插值,interpn,VI=interpn(X1,X2,X3,V,Y1,Y2,Y3,method),多维插值,方法名,说明,nearest,最近插值,linear,线性插值,spline,样条曲线插值,cubic,立方插值,插值函数,Method,选项,34/66,例,7-7,函数,z,=,e,x,+,sin,(,y,)+,y,-1,生成表,7-6,中离散数据，应用不一样插值方法计算在,x,=0.36,，,y,=1.9,处,z,值，并与真值作比较,y,X,0.1,0.2,0.3,0.4,0.5,0.6,0.5,1.0846,1.,1.3293,1.4713,1.6281,1.8015,1.0,1.9466,2.0629,2.1913,2.3333,2.4902,2.6636,1.5,2.6027,2.7189,2.8474,2.9893,3.1462,3.3196,2.0,3.0145,3.1307,3.2592,3.4011,3.5580,3.7314,2.5,3.2036,3.3199,3.4483,3.5903,3.7472,3.92063.0,3.2463,3.3625,3.4910,3.6329,3.7898,3.9632,35/66,例,7-7,插值结果,插值方法,计算结果,真值,3.2796,线性插值,3.2620,三次样条插值,3.2797,最近插值,3.4011,立方插值,3.2784,视频演示,最靠近真实值,36/66,4,非线性方程求解,37/66,4.1,非线性,方程数值求解,常见非线性方程（组）有两种形式,或,直接迭代法、韦斯顿迭代法求解,牛顿迭代法求解,38/66,4.1.1,直接迭代法,先设定,x,初值,x,=,x,0,，代入方程计算,x,1,=,f,(,x,0,),，再把,x,1,作为新初值代入方程计算,x,2,=,f,(,x,1,),直至求得,x,n+1,与,x,n,足够靠近（称为,收敛,），,x,n,即为方程根,直接迭代法求解可写为以下形式,直接迭代法优点是形式简单，易于编程实现。缺点是计算量大、收敛速度慢。普通可经过改进初值、降低收敛要求等方法提升其收敛速度。也可采取其它方法进行求解,39/66,收敛判断准则,绝对偏差,相对偏差,半相对偏差,是用户指定一个很小正数，确定适当,取值有一定难度,优点在于,选取不受方程根数值大小影响。普通取,=0.001,是判断收敛一个很好方法，,取值普通为,0.0001,0.001,40/66,4.1.2,韦格斯顿迭代法,韦斯顿法对直接迭代做了改进，使用前两个计算点信息进行求解。其迭代形式为,韦格斯顿法迭代时需要前面两个计算点数据，可先执行一次直接迭代法计算取得,41/66,4.1.3,牛顿迭代法,又称切线法,若将,f,(,x,),在其根附近进行泰勒级数展开，并取级数线性部分作为,f,(,x,),近似值，可得,由上式可得牛顿迭代公式,如函数一阶导数难以求得，可用差商作为近似导数,42/66,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,函数，给定适当初值对问题求解至关主要，若初值选择不妥，将无法得到正确解。普通可依据经验或简化计算取得适当初值,43/66,例,7-8,试用维里方程计算,200,，,1.013MPa,异丙醇蒸汽摩尔体积,V,与压缩因子,Z,。已知异丙醇维里系数试验值,B=-388cm,3,mol,-1,，,C=-26000 cm,6,mol,-2,视频演示,44/66,例,7-9,600K,下由,CH,3,Cl,和,H,2,O,反应生成,CH,3,OH,，存在以下平衡,CH,3,Cl(g)+H,2,O(g)=CH,3,OH(g)+HCl(g)(1),2CH,3,OH(g)=(CH,3,),2,O(g)+H,2,O(g)(2),已知该温度下,K,p(1),=0.00154,K,p(2),=10.6,。今以等摩尔,CH,3,Cl(g),和,H,2,O(g),开始反应，求,CH,3,Cl,平衡转化率,视频演示,45/66,5,常微分方程（组）求解,46/66,5.1,化工中,常微分方程,(,组,),微分方程中只有一个自变量方程称为常微分方程，自变量个数为两个或两个以上微分方程称为偏微分方程,常微分方程,初值问题,：,给定微分方程及初值条件,边值问题,：,给定微分方程及边界条件,常微分方程（组）解法,解析法和数值法,（惯用）,47/66,初值问题,记为,或,48/66,5.2,常微分方程（组）数值解法,欧拉公式,梯形公式,龙格,-,库塔法,常微分方程组数值解法,49/66,5.2.1,欧拉公式,若常微分方程初值问题求解区间为,，将其等分为,m,步，步长,。记,，对应,x,n,处函数值为,y,n,，则,y,n,可由下式计算,向前欧拉公式,向后欧拉公式,中心欧拉公式,若,y,n+1,同时出现在等号两侧，称为隐式欧拉公式，无法直接求解，普通需采取迭代法计算,50/66,5.2.2,梯形公式,梯形公式也是隐式格式，需要迭代求解,先用显式公式算出初值，再用隐式公式进行一次或数次修正。这一过程称为,预估,-,校正过程,公式为,合并为,51/66,5.2.3,龙格,-,库塔法,工程应用中求解常微分方程最惯用一个有效方法，其计算精度和运算速度较快，易于编程。惯用有二阶、三阶、四阶龙格,-,库塔公式,二,阶龙格,-,库塔公式,，常见形式,或,52/66,三,阶龙格,-,库塔公式,常见形式,53/66,四阶龙格,-,库塔公式,常见形式,54/66,5.2.4,常微分方程组数值解法,将由,m,个一阶方程组成常微分初值问题,其中,可由前边所述解常微分方程各个方法求解,写为向量形式,55/66,5.2.5,高阶常微分方程数值解法,可把高阶常微分方程转化为一阶常微分方程组求解。比如三阶常微分方程,令,将三阶方程,化为一阶方程,56/66,5.3,用,Matlab,求解常微分方程,函数名,求解问题类型,算法,说明,ode45,非刚性问题,Runge-Kutta,一步算法；4，5阶Runge-Kutta方程；累计截断误差达(x)5，精度高，为大部分场所首选算法。,ode23,非刚性问题,Runge-Kutta,一步算法；2，3阶Runge-Kutta方程；累计截断误差达(x)3，计算速度较快，适适用于对精度要求不高情形。,ode23s,刚性问题,Rosenbrock,一步算法；,2,阶,Rosebrock,算法；精度低，若,ode45,失效时，可尝试使用。,ode23t,适度刚性问题,Trapezoidal rule,采取梯形算法求解适度刚性问题。,ode23tb,刚性问题,TR-BDF2,梯形算法，低精度。当精度较低时，计算时间比,ode15s,短,ode15s,刚性问题,NDFs(BDFs),多步算法；,Gears,反向数值微分；精度中等。若,ode45,失效时，可尝试使用,ode113,非刚性问题,Adams,多步,Adams,算法；精度可达,10,-3,-10,-6,。计算时间比,ode45,短,57/66,例,7-10,在间歇反应器中进行液相反应制备产物,B,，其反应网络如图,7-7,所表示。反应温度为,224.6,，反应物,X,大量过剩。各反应均为一级动力学关,系：,各步反应,k,0i,、,E,ai,见表,，试给,出,0-10000,秒各产物浓度改变规律。初始条件为：,t=0,，,CA=1kmol/m3,，,CB=CC=CD=CE=0 kmol/m3,反应网络图,58/66,例,7-10,条件图,参数,取值,参数,取值,k,01,5.7805210,10,E,a1,124670,k,02,3.9231710,12,E,a2,150386,k,03,1.6425410,4,E,a3,77954,k,04,6.26410,8,E,a4,111528,k,05,2.166710,-4,E,a5,0,反应网络图,参数取值,59/66,例,7-10,分析步骤,分析步骤,建立数学模型,建立,odefun,函数,建立主程序,视频演示,60/66,作业,61/66,本章习题,1.,编写程序，实现以下功效：对给定方程组，判断其解情况，若能求解，则求出其解,2.,在学习了怎样判断向量组线性相关之后，怎样判断一个向量能否由其它几个向量线性表示？试编写函数来完成这一功效,62/66,3,3.,判断以下方程组解情况，,若有,解,，求其解,63/66,4.,试依据本章表,7-1,数据应用插值法预测水在,24,、,35,、,48,和,72,时黏度,5.,求非线性方程,在区间,0,1,上根,6.,求非线性方程组,解，初值可设为为,64/66,7,7.,求解以下常微分方程组,65/66,The End,66/66,