Matlab与化学化工计算

《Matlab与化学化工计算》由会员分享，可在线阅读，更多相关《Matlab与化学化工计算（67页珍藏版）》请在金锄头文库上搜索。

1、MatlabMatlab与化学化工计算与化学化工计算本节要点本节要点n n本章背景本章背景n nMatlab基础基础n n方程组求解方程组求解n n数据插值数据插值n n作业作业问题的提出问题的提出n nMATLAB语言与其它语言的关系仿佛和语言与其它语言的关系仿佛和C语言与汇编语言的关系一样语言与汇编语言的关系一样n n计算机语言的发展计算机语言的发展n n标志着计算机语言向标志着计算机语言向“智能化智能化”方向发方向发展，被称为第四代编程语言展，被称为第四代编程语言数值运算数值运算数值运算数值运算解析运算解析运算解析运算解析运算管理、可视化管理、可视化管理、可视化管理、可视化智能化智能化智

2、能化智能化1Matlab基础知识基础知识1.1Matlab简介简介n n1967年由年由ClereMaler用用FORTRAN语言语言设计和编写设计和编写n n1984年年Mathworks公司用公司用C语言完成了语言完成了Matlab的商业化版本并推向市场的商业化版本并推向市场n n经过经过20余年的改进余年的改进，Matlab已发展成为已发展成为一个具有极高通用性的、带有众多实用一个具有极高通用性的、带有众多实用工具的运算平台，成为国际上广泛认可工具的运算平台，成为国际上广泛认可的优秀科学计算软件的优秀科学计算软件Matlab的发展的发展 1984198419841984年，年，年，年，M

3、ATLABMATLABMATLABMATLAB第第第第1 1 1 1版版版版(DOS(DOS(DOS(DOS版版版版) ) ) ) 1992 1992 1992 1992年，年，年，年，MATLAB 4.0MATLAB 4.0MATLAB 4.0MATLAB 4.0版版版版 1994199419941994年，年，年，年，MATLAB 4.2MATLAB 4.2MATLAB 4.2MATLAB 4.2版版版版 1997199719971997年，年，年，年，MATLAB 5.0MATLAB 5.0MATLAB 5.0MATLAB 5.0版版版版 1999199919991999年，年，年，年，

4、MATLAB 5.3MATLAB 5.3MATLAB 5.3MATLAB 5.3版版版版 2000200020002000年，年，年，年，MATLAB 6.0MATLAB 6.0MATLAB 6.0MATLAB 6.0版版版版 2001200120012001年，年，年，年，MATLAB 6.1MATLAB 6.1MATLAB 6.1MATLAB 6.1版版版版 2002200220022002年，年，年，年，MATLAB 6.5MATLAB 6.5MATLAB 6.5MATLAB 6.5版版版版 2004200420042004年，年，年，年，MATLAB 7.0MATLAB 7.0MATL

5、AB 7.0MATLAB 7.0版版版版告告别别DOS版版1993年年MathWorks公司从加拿大公司从加拿大滑滑铁卢铁卢大学大学购购得得Maple的使的使用用权权，推出了，推出了符号符号计计算工具算工具包包5.0的的MATLAB拥拥有有更丰富的数据更丰富的数据类类型和型和结结构、更友构、更友善的面向善的面向对对象、象、更加快速精良的更加快速精良的图图形可形可视视、更广、更广博的数学和数据博的数学和数据分析分析资资源、更多源、更多的的应应用开用开发发工具工具Matlab的优点的优点n n语法简单易学，编程效率高语法简单易学，编程效率高n n高质量、高可靠的数值计算能力高质量、高可靠的数值计算

6、能力n n强大的矩阵运算能力强大的矩阵运算能力n n高级图形和数据可视化处理能力高级图形和数据可视化处理能力n n提供提供600多个常用算法内建函数，以及众多个常用算法内建函数，以及众多面向应用的工具箱多面向应用的工具箱Matlab二维作图二维作图Matlab三维作图三维作图1.2Matlab的界面的界面1.3Matlab的帮助功能的帮助功能n n联机帮助系统联机帮助系统n n命令窗口查询命令窗口查询n nhelphelpn nlookforlookforn n联机演示系统联机演示系统n nDemosDemos“Help”下拉菜下拉菜单中单中“Full Product Family Help”

7、命令打开联命令打开联机帮助系统机帮助系统若不知函数若不知函数确切名，可确切名，可“Lookfor关关键词键词”可查可查helpHelp全部主题全部主题Help指定函数指定函数例例7-1n n查找包含查找包含“diff”关键词的函数关键词的函数n nlookfordifflookfordiffn nSETDIFFSetdifference.SETDIFFSetdifference.n nDIFFDifferenceandapproximatederivative.DIFFDifferenceandapproximatederivative.n nPOLYDERDifferentiatepolyn

8、omial.POLYDERDifferentiatepolynomial.n nDDE23Solvedelaydifferentialequations(DDEs)withDDE23Solvedelaydifferentialequations(DDEs)withconstantdelays.constantdelays.n nDDESDSolvedelaydifferentialequations(DDEs)DDESDSolvedelaydifferentialequations(DDEs)withgeneraldelays.withgeneraldelays.n nDEVALEvaluat

9、ethesolutionofadifferentialDEVALEvaluatethesolutionofadifferentialequationproblem.equationproblem.n n用户输入用户输入的命令的命令查询结果查询结果2 线性方程组求解线性方程组求解2.1线性方程组的一般形式线性方程组的一般形式n n在应用中，常常把线性方程组在应用中，常常把线性方程组n n写成写成AX=b的一般形式，其中的一般形式，其中2.2线性方程组解的判断线性方程组解的判断n n齐次线性方程组齐次线性方程组AX=0，其解的情况可以，其解的情况可以通过系数矩阵通过系数矩阵A的秩和未知数个数的秩和

10、未知数个数n的关的关系来判断系来判断n n如果系数矩阵的秩为如果系数矩阵的秩为如果系数矩阵的秩为如果系数矩阵的秩为n n，方程组只有零解，方程组只有零解，方程组只有零解，方程组只有零解，x x=0=0n n如果系数矩阵的秩小于如果系数矩阵的秩小于如果系数矩阵的秩小于如果系数矩阵的秩小于n n，方程组有无穷多，方程组有无穷多，方程组有无穷多，方程组有无穷多解解解解n n如果系数矩阵的秩大于如果系数矩阵的秩大于如果系数矩阵的秩大于如果系数矩阵的秩大于n n，方程组无解，方程组无解，方程组无解，方程组无解非其次线性方程组解的情况非其次线性方程组解的情况n n非齐次线性方程组非齐次线性方程组AX=b，

11、根据系数矩阵，根据系数矩阵A的秩、增广矩阵的秩、增广矩阵B=Ab的秩和未知数的秩和未知数个数个数n的关系来判断其解的情况的关系来判断其解的情况n n如果系数矩阵如果系数矩阵如果系数矩阵如果系数矩阵A A的秩等于增广矩阵的秩等于增广矩阵的秩等于增广矩阵的秩等于增广矩阵B B的秩且的秩且的秩且的秩且等于等于等于等于n n ，方程组有唯一解，方程组有唯一解，方程组有唯一解，方程组有唯一解n n如果系数矩阵如果系数矩阵如果系数矩阵如果系数矩阵A A的秩等于增广矩阵的秩等于增广矩阵的秩等于增广矩阵的秩等于增广矩阵B B的秩且的秩且的秩且的秩且小于小于小于小于n n ，方程组有无穷多解，方程组有无穷多解，

12、方程组有无穷多解，方程组有无穷多解n n如果系数矩阵如果系数矩阵如果系数矩阵如果系数矩阵A A的秩小于增广矩阵的秩小于增广矩阵的秩小于增广矩阵的秩小于增广矩阵B B的秩，的秩，的秩，的秩，方程组无解方程组无解方程组无解方程组无解例例7-2 判断方程解的情况判断方程解的情况n n解：在解：在解：在解：在MatlabMatlab中输入中输入中输入中输入a=-1-24;211;11-1;a=-1-24;211;11-1;rank(a)rank(a)ans=ans=22齐次线性方程组系数矩阵齐次线性方程组系数矩阵齐次线性方程组系数矩阵齐次线性方程组系数矩阵A A的秩的秩的秩的秩为为为为2 2，小于未知

13、数个数，小于未知数个数，小于未知数个数，小于未知数个数3 3，方，方，方，方程组有无穷多解程组有无穷多解程组有无穷多解程组有无穷多解计算系数矩阵计算系数矩阵A的秩的秩;不能少不能少例例7-2（2）解：解：解：解：a=7028;0281;280196;a=7028;0281;280196;b=1-39-7;b=1-39-7; %b%b为列向量，故输入行向量后转置为列向量，故输入行向量后转置为列向量，故输入行向量后转置为列向量，故输入行向量后转置rank(a)rank(a)%计算系数矩阵计算系数矩阵计算系数矩阵计算系数矩阵A A的秩的秩的秩的秩ans=ans=33rank(ab)rank(ab)%

14、计算增广矩阵计算增广矩阵计算增广矩阵计算增广矩阵AbAb的秩的秩的秩的秩ans=ans=33非齐次线性方程组系数矩阵非齐次线性方程组系数矩阵非齐次线性方程组系数矩阵非齐次线性方程组系数矩阵A A的秩为的秩为的秩为的秩为3 3，增广矩阵的秩为，增广矩阵的秩为，增广矩阵的秩为，增广矩阵的秩为3 3，等于未知，等于未知，等于未知，等于未知数个数数个数数个数数个数3 3，方程组有唯一解。，方程组有唯一解。，方程组有唯一解。，方程组有唯一解。“%”是是Matlab的注释的注释符，符，%后的语后的语句作为注释处句作为注释处理理2.3线性方程组直接求解线性方程组直接求解n n例例例例7-37-3求求求求以下

15、方程组的解以下方程组的解以下方程组的解以下方程组的解n n步骤步骤步骤步骤b b1 1 矩阵除法，验证解矩阵除法，验证解矩阵除法，验证解矩阵除法，验证解n na a判断解的情况判断解的情况判断解的情况判断解的情况b b2 2 逆矩阵法逆矩阵法逆矩阵法逆矩阵法bb3 3rrefrref例例7-4n n求下列方程组的解求下列方程组的解视频视频演示演示3 数据插值数据插值3.1数据插值简介数据插值简介n n在工程领域，许多实验数据常以列表函数或表在工程领域，许多实验数据常以列表函数或表在工程领域，许多实验数据常以列表函数或表在工程领域，许多实验数据常以列表函数或表格的形式存在格的形式存在格的形式存在

16、格的形式存在，如，如，如，如水黏度随温度的列表函数水黏度随温度的列表函数水黏度随温度的列表函数水黏度随温度的列表函数n n在实际使用时，有时需要获得介于表中两个温在实际使用时，有时需要获得介于表中两个温在实际使用时，有时需要获得介于表中两个温在实际使用时，有时需要获得介于表中两个温度结点之间（如度结点之间（如度结点之间（如度结点之间（如1515，2525）的黏度值。而这）的黏度值。而这）的黏度值。而这）的黏度值。而这些数据未在表中出现，需要我们根据已知的数些数据未在表中出现，需要我们根据已知的数些数据未在表中出现，需要我们根据已知的数些数据未在表中出现，需要我们根据已知的数据估算出表中未出现的

17、温度点的黏度数值，这据估算出表中未出现的温度点的黏度数值，这据估算出表中未出现的温度点的黏度数值，这据估算出表中未出现的温度点的黏度数值，这一技术称为一技术称为一技术称为一技术称为插值技术插值技术插值技术插值技术t/0102030405060/mPas1.7881.3051.0040.80120.65320.54920.4698插值的数学定义插值的数学定义n n已知由已知由已知由已知由g g( (X X)()(可能未知或非常复杂可能未知或非常复杂可能未知或非常复杂可能未知或非常复杂) )产生的产生的产生的产生的n n+1+1个个个个离散数据离散数据离散数据离散数据( (x xi i, ,y y

18、i i),i=0,1,2,),i=0,1,2,n n，且，且，且，且这这这这n n+1+1个互异个互异个互异个互异插值结点满足插值结点满足插值结点满足插值结点满足a a= =x x0 0 x x1 1 x x2 2x xn n= =b b，在插值区，在插值区，在插值区，在插值区间间间间 a a, ,b b 内寻找一个相对简单的函数内寻找一个相对简单的函数内寻找一个相对简单的函数内寻找一个相对简单的函数f f( (x x) )，使其，使其，使其，使其满足插值条件满足插值条件满足插值条件满足插值条件f f( (x xi i)=)=y yi i,i=0,1,2,i=0,1,2,n n。再利用已求。再

19、利用已求。再利用已求。再利用已求得的得的得的得的 f f( (x x) )计算任一非插值结点计算任一非插值结点计算任一非插值结点计算任一非插值结点x x* *处的近似值处的近似值处的近似值处的近似值y y* *= =f f( (x x* *) )。其中。其中。其中。其中f f( (x x) )称为插值函数，称为插值函数，称为插值函数，称为插值函数，g g( (x x) )称为被插值称为被插值称为被插值称为被插值函数函数函数函数n n从计算的观点看，插值从计算的观点看，插值从计算的观点看，插值从计算的观点看，插值就就就就是用一个简单函数在是用一个简单函数在是用一个简单函数在是用一个简单函数在某种

20、误差范围内近似的代替原目标函数关系式某种误差范围内近似的代替原目标函数关系式某种误差范围内近似的代替原目标函数关系式某种误差范围内近似的代替原目标函数关系式3.2插值方法插值方法n n线性插值线性插值n n二次插值二次插值n n其他插值方法其他插值方法n n最近（最近（最近（最近（nearestnearest）插值法）插值法）插值法）插值法n n样条曲线（样条曲线（样条曲线（样条曲线（splinespline）法）法）法）法n n埃尔米特（埃尔米特（埃尔米特（埃尔米特（HermiteHermite）法）法）法）法3.2.1线性插值线性插值n n又称两点插值又称两点插值又称两点插值又称两点插值n

21、 n已知两个数据点已知两个数据点已知两个数据点已知两个数据点 x x0 0, ,y y0 0, x x1 1, ,y y1 1(x x0 0 x x1 1) )，求对应，求对应，求对应，求对应于于于于x x( (x x0 0 x x x x1 1) )的的的的y y值值值值解法：由解法：由解法：由解法：由x x0 0, ,y y0 0, ,x x1 1, ,y y1 1构造直构造直构造直构造直线方程线方程线方程线方程并求取在该点的函数值并求取在该点的函数值并求取在该点的函数值并求取在该点的函数值线性插值的优点线性插值的优点是简单，快捷，是简单，快捷，特别是对于插值结点间距较小的特别是对于插值结

22、点间距较小的情况可以取得令人满意的精度情况可以取得令人满意的精度3.2.2二次插值二次插值n n又称拉格朗日三点差值又称拉格朗日三点差值又称拉格朗日三点差值又称拉格朗日三点差值n n根据三个已知点根据三个已知点根据三个已知点根据三个已知点 x x0 0, ,y y0 0, x x1 1, ,y y1 1,x x2 2, ,y y2 2(x x0 0 x x1 1 x x2 2) )，构造二次多项式插值函数构造二次多项式插值函数构造二次多项式插值函数构造二次多项式插值函数y y= =a a0 0+ +a a1 1x x+ +a a2 2x x2 2，并用该函，并用该函，并用该函，并用该函数计算在

23、数计算在数计算在数计算在x x处的处的处的处的y y值值值值n n二次插值公式二次插值公式二次插值公式二次插值公式3.3.1使用使用Matlab进行数据插值进行数据插值n n一维插值一维插值只有一个自变量的插值只有一个自变量的插值n nMatlab提供的一维插值函数是提供的一维插值函数是interp1n n常用语法：常用语法：常用语法：常用语法：YI = interp1YI = interp1( (X,Y,XI,methodX,Y,XI,method) )n n式中式中式中式中n nX X, ,Y Y为已知数据点的为已知数据点的为已知数据点的为已知数据点的x x, ,y y值；值；值；值；n

24、nXIXI 为待插值数据点的为待插值数据点的为待插值数据点的为待插值数据点的x x值；值；值；值；n nYIYI 为返回的插值结果；为返回的插值结果；为返回的插值结果；为返回的插值结果；n n methodmethod用于指定所采用的插值方法用于指定所采用的插值方法用于指定所采用的插值方法用于指定所采用的插值方法插值函数插值函数interp1提供的插值方法提供的插值方法项目项目项目项目取值取值取值取值nearestnearestnearestnearest最近插值最近插值最近插值最近插值linearlinearlinearlinear线性插值，为线性插值，为线性插值，为线性插值，为interp

25、1interp1interp1interp1的缺省方法的缺省方法的缺省方法的缺省方法splinesplinesplinespline分段三次样条插值分段三次样条插值分段三次样条插值分段三次样条插值pchippchippchippchip分段三次分段三次分段三次分段三次HermiteHermiteHermiteHermite插值插值插值插值cubiccubiccubiccubic同同同同pchippchippchippchip不同方法插值的结果不同方法插值的结果n n在在0,2区间内生成区间内生成11个等距的离散点，个等距的离散点，计算函数计算函数y=sin(x)的数值的数值分段三次分段三次分段

26、三次分段三次HermiteHermite插值插值插值插值分段三次分段三次分段三次分段三次样条插值样条插值样条插值样条插值线性插值线性插值线性插值线性插值最近插值最近插值最近插值最近插值例例7-5n n用函数用函数用函数用函数y y= =e ex x生成以下离散数据，使用生成以下离散数据，使用生成以下离散数据，使用生成以下离散数据，使用MatlabMatlab的的的的不同插值方法计算不同插值方法计算不同插值方法计算不同插值方法计算x=2.552.632.772.86x=2.552.632.772.86处的处的处的处的函数值，并与真实值进行比较函数值，并与真实值进行比较函数值，并与真实值进行比较函

27、数值，并与真实值进行比较x2.52.62.72.82.9y12.182513.463714.879716.444618.1741视频视频演示演示例例7-5计算结果的比较计算结果的比较插值方法x2.552.632.772.86真实值12.807113.873815.958617.4615最近插值13.4637 13.4637 16.4446 18.1741 线性插值12.823113.8885 15.9752 17.4823 三次样条插值12.807113.873815.958617.4616分段三次Hermite插值12.806713.873715.958817.4622最接近真实值最接近真实

28、值例例7-6n n已知水在已知水在已知水在已知水在20,21,22,2320,21,22,23的饱和蒸汽压分别为的饱和蒸汽压分别为的饱和蒸汽压分别为的饱和蒸汽压分别为17.54,18.65,19.83,21.07mmHg,17.54,18.65,19.83,21.07mmHg,求求求求20.520.5，21.521.5，22.522.5和和和和2424时水的饱和蒸汽压各是多少？已知时水的饱和蒸汽压各是多少？已知时水的饱和蒸汽压各是多少？已知时水的饱和蒸汽压各是多少？已知2424时水的饱和蒸汽压为时水的饱和蒸汽压为时水的饱和蒸汽压为时水的饱和蒸汽压为22.38mmHg22.38mmHg视频视频演

29、示演示3.3.2多维插值多维插值n n具有多个自变量的插值具有多个自变量的插值n n二维插值二维插值二维插值二维插值n n三维插值三维插值三维插值三维插值n n高维插值高维插值高维插值高维插值函数名典型应用说明interp2ZI = interp2(X,Y,Z,XI,YI,method)二维插值interp3VI = interp3(X,Y,Z,V,XI,YI,ZI,method)三维插值interpnVI = interpn(X1,X2,X3,V,Y1,Y2,Y3,method)多维插值方法名说明nearest最近插值linear线性插值spline样条曲线插值cubic立方插值插值函数插值

30、函数Method选项选项例例7-7n n函数函数函数函数z z= =e ex x+ +sinsin( (y y)+)+y y-1-1生成表生成表生成表生成表7-67-6中的离散数据，中的离散数据，中的离散数据，中的离散数据，应用不同插值方法计算在应用不同插值方法计算在应用不同插值方法计算在应用不同插值方法计算在x x=0.36=0.36，y y=1.9=1.9处的处的处的处的z z值，并与真值作比较值，并与真值作比较值，并与真值作比较值，并与真值作比较yX0.1 0.2 0.3 0.4 0.5 0.6 0.5 1.0846 1.2008 1.3293 1.4713 1.6281 1.8015

31、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.9206 3.0 3.2463 3.3625 3.4910 3.6329 3.7898 3.9632 例例7-7插值结果插值结果插值方法计算结果真值3.2796线性插值3.2620三次样条插值3.2797最近插值3.4011立方插值3.278

32、4视频视频演示演示最接近真实值最接近真实值4 非线性方程的求解非线性方程的求解4.1非线性非线性方程数值求解方程数值求解n n常见的非线性方程（组）有两种形式常见的非线性方程（组）有两种形式或或直接迭代法、直接迭代法、直接迭代法、直接迭代法、韦斯顿迭代韦斯顿迭代韦斯顿迭代韦斯顿迭代法求解法求解法求解法求解牛顿迭代法牛顿迭代法牛顿迭代法牛顿迭代法求解求解求解求解4.1.1直接迭代法直接迭代法n n先设定先设定先设定先设定x x的初值的初值的初值的初值x x= =x x0 0，代入方程计算，代入方程计算，代入方程计算，代入方程计算x x1 1= = f f( (x x0 0) )，再把再把再把再把

33、x x1 1作为新的初值代入方程计算作为新的初值代入方程计算作为新的初值代入方程计算作为新的初值代入方程计算x x2 2= = f f( (x x1 1) )直直直直至求得的至求得的至求得的至求得的x xn+1n+1与与与与x xn n足够接近（称为足够接近（称为足够接近（称为足够接近（称为收敛收敛收敛收敛），），），），x xn n即即即即为方程的根为方程的根为方程的根为方程的根n n直接迭代法求解可写为如下形式直接迭代法求解可写为如下形式直接迭代法求解可写为如下形式直接迭代法求解可写为如下形式直接迭代法的优点是形式简单，易于编程实现。缺直接迭代法的优点是形式简单，易于编程实现。缺直接迭代法

34、的优点是形式简单，易于编程实现。缺直接迭代法的优点是形式简单，易于编程实现。缺点是计算量大、收敛速度慢。一般可通过改进初值、点是计算量大、收敛速度慢。一般可通过改进初值、点是计算量大、收敛速度慢。一般可通过改进初值、点是计算量大、收敛速度慢。一般可通过改进初值、降低收敛要求等方法提高其收敛速度。也可采用其降低收敛要求等方法提高其收敛速度。也可采用其降低收敛要求等方法提高其收敛速度。也可采用其降低收敛要求等方法提高其收敛速度。也可采用其它方法进行求解它方法进行求解它方法进行求解它方法进行求解收敛判断准则收敛判断准则n n绝对偏差绝对偏差n n相对偏差相对偏差n n半相对偏差半相对偏差 是用户指定

35、的一个很是用户指定的一个很小的正数，确定适当的小的正数，确定适当的取值有一定难度取值有一定难度优点在于优点在于的选取不受的选取不受方程根的数值大小的影方程根的数值大小的影响。一般取响。一般取=0.001是判断收敛的一个较好是判断收敛的一个较好的方法，的方法，的取值一般的取值一般为为0.00010.0014.1.2韦格斯顿迭代法韦格斯顿迭代法n n韦斯顿法对直接迭代做了改进，使用前韦斯顿法对直接迭代做了改进，使用前两个计算点的信息进行求解。其迭代形两个计算点的信息进行求解。其迭代形式为式为n n韦格斯顿法迭代时需要前面两个计算点韦格斯顿法迭代时需要前面两个计算点的数据，可先执行一次直接迭代法计算

36、的数据，可先执行一次直接迭代法计算获得获得4.1.3牛顿迭代法牛顿迭代法n n又称切线法又称切线法又称切线法又称切线法n n若将若将若将若将f f( (x x) )在其根附近进行泰勒级数展开，并取在其根附近进行泰勒级数展开，并取在其根附近进行泰勒级数展开，并取在其根附近进行泰勒级数展开，并取级数的线性部分作为级数的线性部分作为级数的线性部分作为级数的线性部分作为f f( (x x) )的近似值，可得的近似值，可得的近似值，可得的近似值，可得n n由上式可得牛顿迭代公式由上式可得牛顿迭代公式由上式可得牛顿迭代公式由上式可得牛顿迭代公式n n如函数的一阶导数难以求得，可用差商作为近如函数的一阶导数

37、难以求得，可用差商作为近如函数的一阶导数难以求得，可用差商作为近如函数的一阶导数难以求得，可用差商作为近似导数似导数似导数似导数4.2用用Matlab求解非线性方程求解非线性方程n nMatlab求解非线性方程求解非线性方程（组）（组）的函数的函数n nfzerofzero：一元非线性方程的求解一元非线性方程的求解一元非线性方程的求解一元非线性方程的求解n nfsolvefsolve：用于非线性方程组的求解用于非线性方程组的求解用于非线性方程组的求解用于非线性方程组的求解n nfzero函数用法函数用法n nx=fzerox=fzero( (funfun, ,x0x0) )n nfsolve函

38、数用法函数用法n nx=fsolve(fun,x0)x=fsolve(fun,x0)fun单变量实值函数，可以是单变量实值函数，可以是Matlab内部函数或用户自定义函数内部函数或用户自定义函数x0若若x0是一个单个的数值，系统是一个单个的数值，系统会将其作为求解的初值，在其附近会将其作为求解的初值，在其附近寻找解；若寻找解；若x0是一个二维向量，且是一个二维向量，且fun(x0(1)和和fun(x0(2)符号相反，符号相反，Matlab将会在将会在x0(1)和和x0(2)区间区间内寻找零点内寻找零点fun用户自定义函数，返回给定变用户自定义函数，返回给定变量量x时方程时方程(组组)的值的值y

39、=fun(x)x0 初值矩阵初值矩阵对于对于fzero和和fsolve函数，给定适当的初值对问题的求解至函数，给定适当的初值对问题的求解至关重要，若初值选择不当，将无法得到正确的解。一般可根关重要，若初值选择不当，将无法得到正确的解。一般可根据经验或简化计算获得合适的初值据经验或简化计算获得合适的初值例例7-8n n试用维里方程计算试用维里方程计算200，1.013MPa的的异丙醇蒸汽的摩尔体积异丙醇蒸汽的摩尔体积V与压缩因子与压缩因子Z。已知异丙醇的维里系数实验值已知异丙醇的维里系数实验值B=-388cm3mol-1，C=-26000cm6mol-2视频视频演示演示例例7-9n n600K6

40、00K下由下由下由下由CHCH3 3ClCl和和和和HH2 2OO反应生成反应生成反应生成反应生成CHCH3 3OHOH，存在，存在，存在，存在下列平衡下列平衡下列平衡下列平衡CHCH3 3Cl(g)+HCl(g)+H2 2O(g)=CHO(g)=CH3 3OH(g)+HCl(g)OH(g)+HCl(g) (1)(1)2CH2CH3 3OH(g)=(CHOH(g)=(CH3 3) )2 2O(g)+HO(g)+H2 2O(g)O(g)(2)(2) 已知该温度下已知该温度下已知该温度下已知该温度下K Kp(1)p(1)=0.00154,=0.00154,K Kp(2)p(2)=10.6=10.6

41、。今以等。今以等。今以等。今以等摩尔的摩尔的摩尔的摩尔的CHCH3 3Cl(g)Cl(g)和和和和HH2 2O(g)O(g)开始反应，求开始反应，求开始反应，求开始反应，求CHCH3 3ClCl的平衡转化率的平衡转化率的平衡转化率的平衡转化率视频视频演示演示5 常微分方程（组）求解常微分方程（组）求解5.1化工中的化工中的常微分方程常微分方程(组组)n n微分方程中只有一个自变量的方程称为微分方程中只有一个自变量的方程称为常微分方程，自变量个数为两个或两个常微分方程，自变量个数为两个或两个以上的微分方程称为偏微分方程以上的微分方程称为偏微分方程n n常微分方程常微分方程n n初值问题初值问题初

42、值问题初值问题：给定微分方程及初值条件给定微分方程及初值条件给定微分方程及初值条件给定微分方程及初值条件n n边值问题边值问题边值问题边值问题：给定微分方程及边界条件给定微分方程及边界条件给定微分方程及边界条件给定微分方程及边界条件n n常微分方程（组）的解法常微分方程（组）的解法n n解析法和数值法解析法和数值法解析法和数值法解析法和数值法（常用）（常用）（常用）（常用）初值问题初值问题n n记为记为或或5.2常微分方程（组）数值解法常微分方程（组）数值解法n n欧拉公式欧拉公式n n梯形公式梯形公式n n龙格龙格-库塔法库塔法n n常微分方程组的数值解法常微分方程组的数值解法5.2.1欧拉

43、公式欧拉公式n n若常微分方程初值问题的求解区间为若常微分方程初值问题的求解区间为若常微分方程初值问题的求解区间为若常微分方程初值问题的求解区间为，将，将，将，将其等分为其等分为其等分为其等分为mm步，步长步，步长步，步长步，步长。记。记。记。记，相应相应相应相应x xn n处的函数值为处的函数值为处的函数值为处的函数值为y yn n , ,，则，则，则，则y yn n可由下式计算可由下式计算可由下式计算可由下式计算n n向前欧拉公式向前欧拉公式向前欧拉公式向前欧拉公式n n向后欧拉公式向后欧拉公式向后欧拉公式向后欧拉公式n n中心欧拉公式中心欧拉公式中心欧拉公式中心欧拉公式若若yn+1同时出

44、现在同时出现在等号的两侧，称为等号的两侧，称为隐式欧拉公式，无隐式欧拉公式，无法直接求解，一般法直接求解，一般需采用迭代法计算需采用迭代法计算5.2.2梯形公式梯形公式n n梯形公式也是隐式格式，需要迭代求解梯形公式也是隐式格式，需要迭代求解梯形公式也是隐式格式，需要迭代求解梯形公式也是隐式格式，需要迭代求解n n先用显式公式算出初值，再用隐式公式进行一先用显式公式算出初值，再用隐式公式进行一先用显式公式算出初值，再用隐式公式进行一先用显式公式算出初值，再用隐式公式进行一次或数次修正。这一过程称为次或数次修正。这一过程称为次或数次修正。这一过程称为次或数次修正。这一过程称为预估预估预估预估-

45、-校正过程校正过程校正过程校正过程n n公式为公式为公式为公式为合并为合并为合并为合并为5.2.3龙格龙格-库塔法库塔法n n工程应用中求解常微分方程最常用的一工程应用中求解常微分方程最常用的一种有效方法，其计算精度和运算速度较种有效方法，其计算精度和运算速度较快，易于编程。常用的有二阶、三阶、快，易于编程。常用的有二阶、三阶、四阶龙格四阶龙格-库塔公式库塔公式n n二二阶龙格阶龙格-库塔公式库塔公式，常见形式，常见形式或或三三阶龙格阶龙格-库塔公式库塔公式n n常见形式常见形式四阶龙格四阶龙格-库塔公式库塔公式n n常见形式常见形式5.2.4常微分方程组的数值解法常微分方程组的数值解法n n

46、将由将由将由将由mm个一阶方程组成的常微分初值问题个一阶方程组成的常微分初值问题个一阶方程组成的常微分初值问题个一阶方程组成的常微分初值问题其中其中可由前边所述的解可由前边所述的解常微分方程的各个常微分方程的各个方法求解方法求解写为向量形式写为向量形式5.2.5高阶常微分方程数值解法高阶常微分方程数值解法n n可把高阶常微分方程转化为一阶常微分方程组可把高阶常微分方程转化为一阶常微分方程组可把高阶常微分方程转化为一阶常微分方程组可把高阶常微分方程转化为一阶常微分方程组求解。例如三阶常微分方程求解。例如三阶常微分方程求解。例如三阶常微分方程求解。例如三阶常微分方程令令令令将三阶方程将三阶方程将三

47、阶方程将三阶方程化为一阶方程化为一阶方程化为一阶方程化为一阶方程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适度刚性问题Trapezoida

48、l rule采用梯形算法求解适度刚性问题。ode23tb刚性问题TR-BDF2梯形算法，低精度。当精度较低时，计算时间比ode15s短ode15s刚性问题NDFs(BDFs)多步算法；Gears反向数值微分；精度中等。若ode45失效时，可尝试使用ode113非刚性问题Adams多步Adams算法；精度可达10-3-10-6。计算时间比ode45短例例7-10n n在间歇反应器中进行液相反应制备产物在间歇反应器中进行液相反应制备产物在间歇反应器中进行液相反应制备产物在间歇反应器中进行液相反应制备产物B B，其，其，其，其反应网络如图反应网络如图反应网络如图反应网络如图7-77-7所示。反应温度

49、为所示。反应温度为所示。反应温度为所示。反应温度为224.6224.6，反应物反应物反应物反应物X X大量过剩。各反应均为一级动力学关大量过剩。各反应均为一级动力学关大量过剩。各反应均为一级动力学关大量过剩。各反应均为一级动力学关系：系：系：系：,各步反应的各步反应的各步反应的各步反应的k k0i0i、E Eaiai见表见表见表见表，试给，试给，试给，试给出出出出0-100000-10000秒各产物的浓度变化规律。初始条件秒各产物的浓度变化规律。初始条件秒各产物的浓度变化规律。初始条件秒各产物的浓度变化规律。初始条件为：为：为：为：t=0t=0，CA=1kmol/m3CA=1kmol/m3，C

50、B=CC=CD=CE=0CB=CC=CD=CE=0kmol/m3kmol/m3反应网络图反应网络图例例7-10条件图条件图参数取值参数取值k015.780521010Ea1124670k023.923171012Ea2150386k031.64254104Ea377954k046.264108Ea4111528k052.166710-4Ea50反应网络图反应网络图参数取值参数取值例例7-10分析步骤分析步骤n n分析步骤分析步骤n n建立数学模型建立数学模型建立数学模型建立数学模型n n建立建立建立建立odefunodefun函数函数函数函数n n建立主程序建立主程序建立主程序建立主程序视频视

51、频演示演示作业作业本章习题本章习题n n1.编写程序，实现如下功能：对给定的编写程序，实现如下功能：对给定的方程组，判断其解的情况，若能求解，方程组，判断其解的情况，若能求解，则求出其解则求出其解n n2.在学习了如何判断向量组线性相关之在学习了如何判断向量组线性相关之后，如何判断一个向量能否由其他几个后，如何判断一个向量能否由其他几个向量线性表示？试编写函数来完成这一向量线性表示？试编写函数来完成这一功能功能3n n3.判断下列方程组解的情况，判断下列方程组解的情况，若有若有解解，求，求其解其解n n4.试根据本章表试根据本章表7-1的数据应用插值法预的数据应用插值法预测水在测水在24、35、48和和72时的黏时的黏度度n n5.求非线性方程求非线性方程在区间在区间0,1上的根上的根n n6.求非线性方程组求非线性方程组的解，初值的解，初值可设为为可设为为7n n7.求解下列常微分方程组求解下列常微分方程组The End