《计算化学》-计算化学-第10章

第10章MATLAB在化学中的应用MATLAB诞生于20世纪70年代，它的编写者是C'leveMoler博士和他的同事。当时，CleveMole:博士与他的同事开发FISPACK和LINPACK的Fortran子程序库，但这两个程序使用存在一定的困难，于是C'leveMole:博士编写了FISPACK和LINPACK的接日程序。CleveMoler给这个接日程序取MATLAB，意为矩阵(MATRIX)和实验室(LABORATORY)的组合，是矩阵实验室的缩写。1982年，CleveMoler和JohnLittle成立了MathWorks公司，正式把MATLAB推向市场，并继续进行MATLAB的开发。下一页返回第10章MATLAB在化学中的应用经过近30年的发展，MATLAB语言已成为当今国际上科学界最具影响力、最有活力的软件，是解决工程与科学实际问题的强大工具。MATLAB擅长于数学计算，可以实现数值分析、优化、偏微分方程数值解的求解问题，而且包含大量的工具箱，可以解决数学和工程领域的绝大多数问题。强大的图形表达功能使该软件不仅可以绘制一般的二维、三维图形，如线条图、条形图、直方图等，还可以绘制工程特性较强的特殊图形，如极坐标图、二维、三维等值线图、三维表面图、流线图等，此外还可以生成快照图和进行动问制作。可扩展功能用户可以根据自己的需要编写文件，组成自己的工具箱，以解决本领域内常见的计算问题。MATLAB的优势主要体现在以下几个特点:①语言简洁紧凑，使用方便灵活，库函数极其丰富，压缩了一切不必要的编程工作;上一页下一页返回第10章MATLAB在化学中的应用②运算符丰富。MATLAB是由C语言编写的，它提供了和C语言几乎一样的运算符;③语法限制不严格，程序设计自由度大;④程序可移植性很好，基本上不用做修改就可在各种型号的计算机和操作系统上运行;⑤图形功能强大，数据的可视化非常简单;⑥功能强劲的工具箱是MATLAB的另一重大特色，常用的工具箱有:分析与综合工具箱(AnalysisandSynthesisToolbox)、神经网络工具箱(NeuralNetworkToolbox)等;⑦原程序的开放性，除内部函数以外，所有的MATLAB的核心文件和工具箱文件都是可读可改的源文件，用户可通过对源文件的修改以及加入自己的文件构成新的工具箱。上一页下一页返回第10章MATLAB在化学中的应用目前，MATLAB可以用来进行:①数值分析;②数值和符号计算;③工程和科学绘图;④控制系统的设计与仿真数字图像处理;⑤数字信号处理;⑥通信系统设计与仿真;⑦财务与金融工程。简言之，MATLAB具有非常强大的数值计算功能、图形可视化功能和符号运算功能，而且简单易学，扩展性好，可以与其他面向对象的高级语言进行混合编程。上一页返回10.1MATLAB编程基础10.1.1MATLAB的基本操作1.MATLAB的帮助功能(1)用help命令在命令窗口中只输入help命令，则可以列出MATLAB扫搜路径的所有目录:>>help如果知道函数名，如corvefit，则可用下列命令了解其用法:>>helpcorvefit如果用目录名代替函数名，则列出目录内容，例如:>>helpmatlab\general当不知道某函数的确切名字时，先用lookfor命令扫搜，然后用help:>>lookforcurv下一页返回10.1MATLAB编程基础(2)联机帮助在MATLAB界面中单击工具条的问号(?)按钮或单击help菜单中中MATLABHelp选项，即可打开联机帮助界面。(3)演示帮助——看范例>>demo(4)PDF帮助文档完整的MATLAB软件带有PDF文件类型的帮助文档，每个MATLAB工具箱都有相应的PDF帮助文档，可以用AdobeAcrobatReader软件打开进行查阅。(5)MATLAB网站资源访问MathWorks公司的网站http://。上一页下一页返回10.1MATLAB编程基础2.MATLAB的基本操作命令(1)操作符算术操作符:+-*/^(乘方)\(左除)关系操作符:=(等于)~=(不等于)逻辑操作符:&(逻辑与)|(逻辑或)~(逻辑否)(2)MATLAB工作环境用who或whos可以观察当前工作环境中的变量:>>who%只查看变量名>>whos%详细了解变量名、大小、所占内存空间、数据类型等用size命令可以了解矩阵的大小:>>a=[123;456]>>[m,n]=size(a)%m代表行数，n代表列数上一页下一页返回10.1MATLAB编程基础用length命令可以了解向量的大小:>>v=[678]>>p=length(v)%这里p代表向量的长度用clear命令可以清除工作环境中的全部变量和函数，以释放相应的内存空间:>>clear，who也可以用clear命令从内存中清除指定的变量:>>clearx，who%从内存中清除变量x用clc命令可以清屏并使光标移动到屏幕最左上方:>>clc上一页下一页返回10.1MATLAB编程基础MATLAB所有计算都是按双精度进行的，但计算结果有多种显示格式，默认格式以5个数字显示，也可以按格式命令(format)指定显示格式:>>a=sqrt(2)>>formatshort，a>>formatlonge，a在工作环境窗口可以执行DOS和UNIX命令，调用格式分别是:Dos('command')和unix(’command')，如:>>dos('c:')，unix('c:')>>dos('dir')，unix('dir')>>dos('cdC:\MATLAB7.0\toolbox')，unix('cdC:\MATLAB7.0\toolbox')上一页下一页返回10.1MATLAB编程基础10.1.2MATLAB的变量、常数和数据类型1.MATLAB的变量及其作用域与其他汇编语言一样，MATLAB采用变量来存储数值。但是，MATLAB不需要声明，可以直接使用。根据变量作用域的不同，可以将变量分为局域变量和全局变量。在默认的情况下，函数内的变量属于局域变量，它只在函数内有效，而在该函数外部是不可用的。全局变量对整个程序的所有函数都有效，全局变量可以用global关键词定义。2.MATLAB的数组、向量和矩阵MATLAB数组是MATLAB语言唯一能处理的对象类型，MATLAB所有类型的变量包括标量、向量、矩阵、字符串、单兀数组、结构体和对象等，都以MATLAB数组方式储存。上一页下一页返回10.1MATLAB编程基础MATLAB数组的储存方式:与FORTRAN语言一样，所有MATLAB数据按下列方式存储。如在命令窗可键入如下矩阵a:>>a=['house'；'floor'；'porch']a=housefloorporch>>size(a)ans=35向量、矩阵的操作:由于矩阵使用最为广泛，MATLAB的基本实体是矩阵。用中括号可以引入一个矩阵，其中，一行中的元索用空格或逗号隔开，而行与行之间用分号或回车分开，如:>>t=[123;256]上一页下一页返回10.1MATLAB编程基础10.1.3数据的输入和输出1.MATLAB数据的输入在MATLAB编程中，有多种方法可以将数据输入到MATLAB程序中。(1)利用M文件产生数据文件利用文本编辑器可产生一个数据M文件，用于保存已知参数。这样，可以在MATLAB程序中调用该数据M文件，将有关变量及其数据直接调入MATLAB内存中。(2)用Load命令从MAT文件或文本文件中读取数据用Load命令可以将指定文件的数据装载到MATLAB工作环境中，若装载文件扩展名为.mat，则装载时将显T以Save命令保存的变量。上一页下一页返回10.1MATLAB编程基础(3)用fscanf函数用fscanf函数可以从键盘或文件读取格式化数据，所用的控制符及用法与函数fprintf()完全类似。(4)用提示输入函数input2.MATLAB数据的输出(1)用Save命令用Save命令可以将MATLAB工作环境中的几个或全部变量保存到文件中。(2)用fprintf函数用fprintf函数可以按格式将数据输出至屏幕或写格式化数据到文件中。该函数与C语言的同名函数功能和用法大致相同。(3)用函数disp()将结果输出至屏幕上一页下一页返回10.1MATLAB编程基础10.1.4图形可视化MATLAB具有很强的图形处理功能，可以绘制二维图形、三维图形和四维表现图等，还能对图形进行着色、消隐、光照、渲染等处理。1.二维图形利用plot和line函数可以在二维平面上绘制二维图形。其中，line是直接对图形进行操作的底层画图命令，plot是在底层命令基础上建立起来的高层绘图命令。上一页下一页返回10.1MATLAB编程基础如在窗口中输入命令:x=0:pi/18:2*pi;%给出0~2π区间，间隔为π/18的横坐标y=sin(x);%计算出纵坐标plot(x,y);%绘制图形(如图10.1所示)2.三维图形三维图形主要包括二维曲线和曲面图，前者用plot3函数绘制，后者用mesh或surf绘制。上一页返回10.2MATLAB在化学、化工中的应用MATLAB的函数库十分庞大，包含2000多个函数，几乎包括了一般用户所需的所有函数。与MATLAB配用的SymbolicMathToolbook2.0允许用户在三个不同层次上进行符号计算。①在进行符号对象定义后，直接利用MATLAB格式进行矩阵分解、微分、积分、积分变换、代数方程求解、微分方程求解等运算;②借助maple指令，把单个MAPLE格式的指令送往MAPLE引擎计算;③借助procread把整段MAPLE程序送往MAPLE计算。利用它们可以解决在化学化工中棘手的复杂数学问题。下一页返回10.2MATLAB在化学、化工中的应用10.2.1数据插值插值在科学研究和工程实践中应用广泛，数据插值是利用一批已知的实测数据，采用某种算法平滑地估计出测量点之间的数据。常见的插值方法有Lagrange多项式插值、牛顿插值、分段插值,Hermite插值和三次样条插值，MATLAB主要提供分段线性插值和三次样条插值，而对于Lagrange多项式插值、牛顿插值、Hermite插值，MATLAB没有提供相应的函数，想用这几种函数可以自己编写函数。上一页下一页返回10.2MATLAB在化学、化工中的应用【例10.1】青霉素发酵时i司数据见表10.1,请用插值方法估算t=10,30,50,70,90,110,130,170,190时青霉索的浓度。程序过程:clearallclc%调入青霉索发酵的实验数据t=[020206080100120120160180200];%发酵时间:c=[010616003000515886009230109501028096209200];%青霉索浓度ti=[10,30,50,70,90,110,130,170,190];%求插值点出青霉索浓度ti=10:20:190ci=spline(t，c，ti)PP=spline(t，c)tj=0:0.5:200;上一页下一页返回10.2MATLAB在化学、化工中的应用cj=ppval(pp，tj);%图形输出plot(t，c，'ko'，ti，ci，'b*'，tj，cJ，'b-')legend('实验值'，'插值点'，'样条插值曲线')xlabel('时间(h)')ylabel('青霉索浓度')计算结果:ti=1030507090110130170190ci=1.0e+0.04*-0.02740.03330.22230.42400.74200.90611.02671.08220.98610.9489第一个值为负值，不符合实际物理意义，这是由于插值误差造成的。观察得出的插值图形，如图10.2所示，结果比较令人满意。上一页下一页返回10.2MATLAB在化学、化工中的应用10.2.2方程(组)求解1.一元非线性方程数学模型:非线性函数f(x)=0。函数fzero调用格式:[x,fval]=fzero(@myfun,x())功能:找出函数值为零的点。说明:x()为初始值。【例10.2】对气液平衡物系，由于组分的饱和蒸气压与温度的关系通常为非线性的，若给定液相或气相组成和温度，一般需要试差才能确定与之相平衡的液相或气相组成和温度。例如常压下苯(A)和甲苯(B)物系，溶液中苯的物质的量为0.2mol/L，利用平衡关系式:上一页下一页返回10.2MATLAB在化学、化工中的应用利用安托因方程:而安托因公式是只含有一个末知数(泡点温度)的非线性方程。将公式变换成等于零的函数式。上一页下一页返回10.2MATLAB在化学、化工中的应用①编写名为myfnn.m的M文件。functionF=myfun(x)F=(106.7-10^(6.080-1345/(x+219.5)))/(10^(6.030-1211/(x+220.8))-10^(6.080-1345/(x+219.5)))-0.20存盘。上一页下一页返回10.2MATLAB在化学、化工中的应用②在命令窗口输入:x0=90;%泡点温度的初始值[x，favl]=fzero(@myfun，x())回车，在命令窗口立即得到:F=0.4415F=0.5433F=1.5821e-015F=-2.7756e-016x=103.9176favl=-2.7756e-016即泡点温度为103.9℃时，函数值fval=-2.7756×10-16，满足要求。上一页下一页返回10.2MATLAB在化学、化工中的应用③在命令窗口输入:P0A=10^(6.030-1211/(x+220.8))；yA=P0A*0.20/106.7回车，并得到:yA=0.3745可见，泡点温度为103.9℃时，对应的气相组成yA=0.3745，省去了试差的烦琐，提高了计算效率。上一页下一页返回10.2MATLAB在化学、化工中的应用2.非线性方程组数学模型:含有n个末知数与n个方程的非线性方程组F(x)=0函数solve调用格式:g=solve(eq1,eq2,...,eqn)函数fsolve调用格式:x=fsolve(fun，x())[x,fval]=fsolve(fun,x())[x,fval,exitflag]=fsolve(fun,x(),option,P1,P2...)【例10.3】细胞反应得率系数计算。在分批培养中，培养基组成时刻在变化，某一时刻的瞬时得率

需计算葡萄糖生成乙醇发酵中的得率系数。该发酵过程限制性培养基中细胞生长方程、产物生成量计量式及底物比消耗速率可分别表示如下(XATP，代表能量的消耗与生成)。上一页下一页返回10.2MATLAB在化学、化工中的应用①②(假定产物除CO2外只有乙醇，且生成的ATP全部用于细胞生长过程的耗能，即不考虑维持代谢)。③该问题的实质是求解非线性方程组。由得率因子方程:式中，S→X反应的化学计量系数VS/X=1.12/1;细胞的摩尔质量为24.2;葡萄糖的摩尔质量为30。上一页下一页返回10.2MATLAB在化学、化工中的应用得由得设则上述问题变为求解下列非线性方程组:调用MATLAB的solve函数，采用g=solve(eql,eq2,...,eqn)的形式。在命令窗口输入:结果:x1=-5.1234,0.3698;x2=4.4922,0.3905则上一页下一页返回10.2MATLAB在化学、化工中的应用10.2.3微分方程(组)求解研究反应动力学的数学本质就是解微分方程或微分方程组的问题。MATLAB数值计算方面的最大变化是增添了一组常微分方程(ODE)数值解算程序。这组解算程序无论是算法还是软件结构都十分精良，它包含ode23,ode45,odell3,ode23t,odel5s,ode23s,ode23tb等不同解算指令，用以解算包括stiff方程在内的各种微分方程。MATLAB为解常微分方程问题所设计的文件十分严整，包括解算指令solver、被solver调用的微分方程描述文件、进行积分算法参数设置的odeset和odeget、解算输出指令odeplot,odephas2等。上一页下一页返回10.2MATLAB在化学、化工中的应用1.微分方程求解调用格式:diff(S,'v')%求符号表达式S对指定的独立变量v的微分

diff(S,n)%参数n为指定求微分的阶数

diff(S,'v',n)【例10.4】研发某新型抗生索农药，通过单因子试验与均匀试验，筛选活性组分最佳(以抑菌圈直径大小衡量)的发酵培养基配方。在试验中，以配方用量为自变量，实测效价为因变量，通过二次多项式回归分析，建立了多元二次回归数学模型:将y分别对x1、x2、x3、x4、x5求偏导，当函数为零时，应有极大值。上一页下一页返回10.2MATLAB在化学、化工中的应用在MATLAB窗口输入:symsx1、x2、x3、x4、x5；f=sym('0.06541*x12-4.948*x22-48.473*x32+1.427*x42-44.658*x52+7.348*x2*x3-0.352*x1+12.371*x2-1.254*x3-0.397*x4+14.171*x5+33.576');A=diff(f，x1)H=diff(f，x2)C=diff(f，x3)D=diff(f，x4)E=diff(f，x5)上一页下一页返回10.2MATLAB在化学、化工中的应用输出为得到5个线性方程组上一页下一页返回10.2MATLAB在化学、化工中的应用在MATLAB窗口输入:输出为方程最佳值依次为上一页下一页返回10.2MATLAB在化学、化工中的应用将最优解代入上述方程，在MATlAB窗口输入:优化的培养基发酵后，抑菌圈直径应为42.27mm。按照优化条件进行试验，其结果为43.25mm，在预测范围内，所得到的优化发酵培养基的最佳配方为:淀粉2.69%，花生饼粉1.31%，(NH4)2SO40.08%，CaCO30.74%，NaCl0.16%。上一页下一页返回10.2MATLAB在化学、化工中的应用2.微分方程组求解解常微分方程组需调用ODE函数:ode23,ode15,odell3,ode23t,odel5s,ode23s,ode23tb等。调用格式:ode45(odefun,bspan,x())[t,x]=ode45(odefun,bspan,x())上一页下一页返回10.2MATLAB在化学、化工中的应用【例10.5】细胞反应动力学过程计算。酵母在葡萄糖上生长并生成乙醇过程动力学可以用下述公式表示:该问题的实质是求解常微分方程组，特征常数及初始条件已给定，可进行数值求解。用MATLAB四~五阶龙格——库塔法函数ode45求解常微分方程组。上一页下一页返回10.2MATLAB在化学、化工中的应用令cx=x(1),cp=x(2),cs=x(3)。在MATLAB环境下编程:①编写名为dlx.m的M文件。Functiondxdt=dlx(t.x)②在命令窗口输入:tspan=[0100];%时间间隔x0=[0.10100];%cx、cp、cs的初始值ode45(@dlx，tspan，x())得到酵母、乙醇和葡萄糖各组分随时间的变化曲线，如图10.3所示。上一页下一页返回10.2MATLAB在化学、化工中的应用③若在命令窗口输入:tspan=[020206080100];%时间间隔x0=[0.10100];%cx、cp、cs的初始值

[t，x]=ode45(@dlx，tspan，x())结果为t=020406080100上一页下一页返回10.2MATLAB在化学、化工中的应用得到发酵0h,20h,10