工程应用数学- 应用实例

工程应用数学数学计算方法与软件的工程应用,本学期主要教学内容,部分应用实例矩阵分析（Chap2）、线性方程组（Chap3）、数值逼近（Chap5）、非线性方程组实验设计与数据处理实验设计与分析（Chap7、8）模型参数回归（Chap8）神经网络（Chap8）其他模拟退火、遗传算法(Chap8)微分方程组（Chap9、10）参数估计,工具软件：MATLABStatistica,数学模型,数学模型是系统某种特征本质的数学表达式即用数学式子（如函数式、代数方程、微分方程、积分方程、差分方程等）来描述（表达、模拟）所研究的客观对象或系统在某一方面的存在规律。用模型描述实体的过程称为建模或模型化。理想的数学模型必须满足以下两点:可靠性：在允许的误差值围内，它能反映出该系统的有关特性的内在联系。适用性：它须易于数学处理和计算。,建立数学模型的一般方法,建立模型的方法大致有两种：实验归纳法理论分析法建立模型的一般步骤：通过对系统的分析，根据问题的性质和精度的要求，作出合理性假设、简化。抽象出系统的物理模型。在此基础上确定输入、输出变量和模型参数，建立数学模型。一般来说，在不降低精度的条件下,模型变量的数目越少越好。检验和修正所得模型。手段是将模型计算结果与实验结果做对比.模型含有无关或关系不大的变量；模型遗漏了重要的有关变量；模型参数不准确；数学模型的结构形式有错；模型反映系统的精确度不够。,数值计算在化学化工中的应用,计算无处不在！,模拟-实验-理论的关系,由模拟所得到的结果，再经科学实验得到的实验值相互比较，最后经过理论的验证，便可以得到一个真正完美的结论。所以在近代科学研究中，必须以实验验证理论，以模拟辅助实验，三者并行，如此便能达到相辅相成的效果。,提高数学建模效率的方法,建立模型,调用算法,分析结果,低效率,高效率,瓶颈,重点,重点,专有的商业软件,实质：数学模型和计算方法的有机集成过程模拟：ASPENPLUS、PRO/II、HYSYS、gPROMS等分子模拟：Gaussian、Cerius2、MaterialsStudio、InsightII、DISCOVERYSTUDIO、SYBYL、HyperChem、CHEMOFFICE等计算流体动力学：CFX、FLUENT、StarCD等优势：技术成熟、系统稳定、资料丰富、技术交流方便。缺陷：价格高，解决对象为已有的成熟的工程问题，缺少新的研究课题的数学模型。,对于科学研究领域，通过建模、编程解决新的模型问题成为必要。,数值计算的工具,程序设计语言BASIC/VisualBASIC(VB)PASCAl/DelphiC/C+(VisualC+、C+Builder)FORTRAN/CompacVisualFortran6.6(CVF)数学软件包Mathematica（数学演算）MathCADMaple（符号计算）MATABSAS、SPSS、STATISTICA（统计分析）,执行效率高、有丰富的数值计算源程序或库文件，如NumericalRecipes、IMSL库以及网络资源NetLib。编程能力要求高。,算法齐全，计算、图形可视化和符号运算功能强大，且简单易学、扩展性好，也支持与其他高级语言混合编程。既是专业数学软件，又是一种编程语言，编程效率高，且代码公开。内建丰富的函数和工具箱。得到第三方公司的开发支持,MATALB在化学化工中的应用,数值计算AlkisConstantinides,navidMostoufi.NumericalMethodsforChemicalEngineeringwithMATLABApplications.PrenticeHall,1999MichaelB.Cutlip102-101000;01-110-1010;0000.5.5000;-1-20-1-100-20;120110-241;0000001-1-2;000000001;B=svd(C),B=6.39463.26802.77861.55790.88390.00000.00000.0000,结果表明矩阵的秩为5，故此体系的独立反应数目为5。,若需要求构成1组独立反应的反应式，如何做？,反应矩阵转置,求解,对所有可能的独立反应组进行全范围搜索，然后判断其秩是否等于反应矩阵的秩，即可求得独立反应的构成。,C=-1-1-1000-100;102-101000;01-110-1010;0000.5.5000;-1-20-1-100-20;120110-241;0000001-1-2;000000001;A=C;a=rank(A);disp(sprintf(n最高秩=%i,a)disp(sprintf(n可形成rank=%i的有下列几组独立反应方程式,a),fori=1:length(A)forj=i:length(A)B=A(i,:);A(j,:);A(k,:);A(l,:);A(m,:);b=rank(B);ifb=adisp(sprintf(%i%i%i,i,j,k,l,m)elseendend,作业,选择一个复杂反应体系，分析有多少个独立的反应数？并求构成1组独立反应的反应式（方程式号码）。,作业格式要求：,问题说明，包括背景、要解决的问题问题的分析，提出数学模型计算机求解的方法、计算流程、源程序结果分析与讨论,Ex2：线性方程组的应用,求解线性方程组的方法：直接解法：Gauss消元法；矩阵分解方法；Cholesky分解函数chol；A=LLTLU分解函数lu；A=LU正交三角分解函数qr等。A=QR迭代解法：Jacobi迭代法Gause-Seidel迭代法SuccessiveOver-Relaxation(SOR)迭代法特殊线性方程组：稀疏线性方程组共轭梯度迭代法；三对角方程组追赶法,MATLAB求解线性方程组方法,用矩阵左除操作符“”求解；x=Ab用函数solve计算解析解；g=solve(eq1,eq2,eqn,var1,var2,varn)其中eq1,eq2,eqn是方程组的n个方程；var1,var2,varn是方程的变量列表。通过矩阵求逆函数inv；x=inv(A)*b矩阵分解方法；Cholesky分解函数chol；L=chol(A);y=Lb;x=LyA=LL;Ly=b;Lx=yLU分解函数lu；L,U,P=lu(A);y=Lp*b;x=UyLU=pA;Ly=pb;Ux=y正交三角分解函数qr等；迭代法和追赶法等须自行编制程序。,乙醇精馏过程的物料平衡,生产工业乙醇的二级精馏过程见下图。设过程处理量10000kg/hr，组成为80水，10醇和10有机物(重量)，精馏塔回流比为3。第一精馏塔顶产物含醇60，而第二精馏塔顶产物含醇95。第一精馏塔底含进料有机物的80，剩余的在第二塔底，两塔塔底物料中都不含乙醇。求解每股物料的量。,线性方程组问题,W表示水A表示乙醇R表示有机物料,包含再沸器和冷凝器的第一蒸馏塔,蒸馏塔,回流比,塔底物料中都不含乙醇,精馏塔回流比为3,按照回流比3计算，W3/W5=4/3,包含再沸器和冷凝器的第二蒸馏塔,蒸馏塔,回流比,塔底物料中都不含乙醇,精馏塔回流比为3,按照回流比3计算，W9/W11=4/3,约束方程,在列方程时要考虑对角元素不为零。较方便的是从第一流股开始，顺序进行，这里从3流股开始，13流股结束，矩阵表示见下图。,此外还有三个独立的线性约束：,第一精馏塔顶产物含醇60,第二精馏塔顶产物含醇95,第一精馏塔底含进料有机物的80,矩阵形式,在列方程时要考虑对角元素不为零。较方便的是从第一流股开始，顺序进行，这里从流股3开始，流股13结束,源程序数据输入,%SetvaluesforAandB%N=19;A=zeros(N,N);B=zeros(N,1);A(1,4)=-4/3;A(2,5)=-4/3;A(3,6)=-4/3;A(4,9)=-3;A(5,10)=-3;A(6,11)=-3;A(7,9)=1;A(9,10)=-2/3;A(9,11)=1;A(11,8)=1;A(12,14)=-4/3;A(13,15)=-4/3;,A(14,18)=-3;A(15,19)=-3;A(16,9)=-1;A(16,18)=1;A(17,11)=-1;A(18,19)=-5/95;A(19,10)=-1;forI=1:NA(I,I)=1;%B(I)=0;end;B(7)=8000;B(8)=800;B(10)=1000;B(11)=1000;,源程序求解与结果,%directsolve%X=AB;%LUfactorL,U,P=lu(A);Y=LP*B;X=UY;%resultdisplaydisp(sprintf(%8.1f,X),ethanol1866.74000.0800.01400.03000.0600.07533.3800.0466.71000.0200.0210.54000.0157.93000.0414.0200.052.61000.0,物料流股数据,Ex3：数值逼近的应用,插值拟合数值积分数值微分,MATLAB的插值函数,一般调用函数interp1、interp2和interp3进行多项式插值计算；yi=interp1(x0,y0,xi,method)其中：x0和y0分别为观察节点和观察值向量；xi表示待求函数值的节点；method指定插值的方法，取值为linear时表示进行线性插值，spline时表示进行三次样条插值，cubic时表示进行三次多项式插值；输出yi表示在点xi处的多项式插值的值。调用函数spline进行多项式样条插值计算；y=spline(x0,y0,x)或pp=spline(x0,y0)其中：x0和y0分别为观察节点和观察值的向量；前者输出在点x处的三次样条值，后者得到pp格式(PiecewisePolynomial)的样条函数形式。fnbrk(pp)输出样条多项式格式,插值方法比较,例：不同插值方法的比较做一个函数sin(x2)在区间0，2上40个函数值的表：用interp1来计算中间点的函数值，而不用sin(x)求：用样条插值来计算中间点的函数值：用sin计算得到的正确结果,x=linspace(0,2*pi,40);y=sin(x.2);,values=00.604980.35591,values=interp1(x,y,0pi/23),better=interp1(x,y,0pi/23,spline),correct=sin(0pi/23.2),better=00.624140.40982,correct=00.624270.41212,样条插值可得更高精度的结果，在表中用更多的数据点也可以使精度更好。,三次样条数据插值：示例,例：对离散地分布在y=exp(x)sin(x)函数曲线上的数据点进行样条插值计算,x=024581212.817.219.920;y=exp(x).*sin(x);xx=0:.25:20;yy=spline(x,y,xx);plot(x,y,o,xx,yy)pp=spline(x,y);fnbrk(pp),Theinputdescribesappformbreaks(1:l+1)coefficients(d*l,k)piecesnumberlorderkdimensiondoftarget,MATLAB的插值函数,专门进行样条插值的样条工具箱(SplineToolbox)；函数csape可以构造各种边界条件下的三次样条插值；函数csapi用于确定not-a-knot边界条件下的三次样条插值；函数csaps可产生不同光滑程度的三次样条插值；函数cscvn则构造周期三次样条插值。如函数csape：pp=csape(x0,y0,conds,valconds)其中conds和valconds用于描述各种边界条件,MATLAB的拟合函数,函数polyfit可进行离散数据的多项式拟合；p=polyfit(x0,y0,n)其中：x0和y0分别为观察节点和观察值向量；n表示插值多项式的次数；输出值p表示插值多项式的系数；函数lsqcurvefit用于处理非线性模型的最小二乘拟合；beta=lsqcurvefit(fun,beta0,xdata,ydata)其中：xdata和ydata分别为观察节点和观察值向量；fun为自定义的非线性拟合模型；beta0为拟合参数的初始值；输出值beta为拟合模型中的参数。,MATLAB的数值积分函数,用梯形积分的函数trapz进行计算；采用自适应方法的积分函数quad（低阶方法，自适应辛普森方法）和quad8（高阶方法，自适应柯特斯公式计算）计算；Q=quad(fun,a,b,tol)或Q=quad8(fun,a,b,tol)其中：参数fun代表被积函数；a和b分别表示定积分的下限和上限；tol则表示绝对误差限，默认值为1e-6；输出值Q就是定积分的近似值。符号积分的函数int；R=int(fun,v)和R=int(fun,v,a,b)其中：fun为积分函数；v表示积分变量；a和b分别表示定积分的下限和上限；输出值根据调用形式为不定积分或定积分值。,积分方法示例,Trapz梯形积分自适应Simpleson积分法,X=-1:.1:1;Y=1./(1+25*X.2);T=trapz(X,Y),T=0.54924,fun=inline(3*x.2./(x.3-2*x.2+3);Q1=quad(fun,0,2)Q2=quadl(fun,0,2),Q1=3.7224Q2=3.7224,数值积分方法比较,（1）将(0，pi/2)10等分，步长为pi/20，利用trapz(x)函数按梯形算法计算（2）使用quad(fun,a,b)函数命令的辛普森公式法：（3）蒙特卡罗均值估计法：蒙特卡罗法是一种随机实验方法。本例使用sum(x)函数，sum(x)函数的功能是输入数组x，输出x的和。,z=0.99794,m=pi/20;x=0:m:pi/2;y=sin(x);z=trapz(y)*m,z=quad(sin,0,pi/2),r=100000;x=rand(1,r);y=sin(x*pi/2);z=sum(y)*pi/2/r,z=1,z=1.001,上面三种方法得出的结果不同这是由于所选方法各自产生的误差不同所造成的。,MATLAB的数值微分函数,函数diff进行计算；函数gradient计算多元实值函数的梯度；函数jacobian计算多元向量值函数的雅可比矩阵。,例：汽液平衡的一致性检验计算,准确的汽液平衡数据可以用来计算液相活度系数和Gibbs过量自由能。对下表的双组分汽液平衡数据，计算液相活度系数，并对结果进行热力学的一致性检验，假定气相是理想的。,数据来源：S.WeissmanandS.E.Wood,Vapor-LiquidEquilibriumofBenzene-2,2,4-trimethylpentaneMixtures,J.Chem.Phys.,vol.32,1960,p.1153,汽相摩尔分率,液相摩尔分率,理论分析,对汽液平衡数据进行热力学一致性检验可用Gibbs-Duhem方程。对双组分情况，Gibbs-Duhem方程的形式是：如采用面积检验法，有：,活度系数计算,无限稀释的活度系数估算(1),无限稀释的活度系数估算(2),无限稀释的活度系数估算(3),无限稀释的活度系数估算(4),一致性检验方法,作出整个区间0,1的待积分图形：,一致性检验：,源程序（1）,插值画图,functionvle_spline%vle_spline.m%clearallclcglobalsp;x0=0.0;x1=1.0;xi=0.00.08190.21920.35840.38310.52560.84780.98721.0;yi=-0.3709-0.3386-0.2840-0.1787-0.1670-0.03670.38760.67670.6901;sp=spline(xi,yi);,figure(1);holdon;AXIS(01-11);%x4plot=linspace(xi(1),xi(end),200);y4plot=fnval(sp,x4plot);plot(xi,yi,o,x4plot,y4plot,-)%line(xp(i-1)xp(i),yp(i-1)yp(i)line(0,1,0,0);xlabel(x1);ylabel(Activityratio);holdoff;,源程序（2）,N=quadl(func1,x0,x1);disp(计算结果:)fprintf(n积分区间为(%f,%f)n,x0,x1)fprintf(积分结果为：%fn,N),计算结果:积分区间为(0.000000,1.000000)积分结果为：0.007272,面积之和还是面积之差？,源程序（3）,sp=spline(xi,yi);y2=abs(yi);sp2=spline(xi,y2);N1=quadl(func1,x0,x1);%积分结果N2=quadl(func2,x0,x1)disp(计算结果:)fprintf(n积分区间为(%f,%f)n,x0,x1)fprintf(面积之差为：%fn面积之和为：%fn,abs(N1),abs(N2)per=abs(N1/N2)*100;fprintf(面积之差与面积之和的百分比为：%fn,per)ifperroots(230101),代数方程的图解法,采用函数ezplot(f)画出表达式f=f(x)的图，默认的范围是：-2x2。也可以指定范围，如ezplot(f,min,max)画出f=f(x)的图，区域范围是：minholdoffezplot(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5,05)holdonline(0,5,0,0)grid,t=3.52;x=exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5,x=3.1545e-004,图解法由于受到空间维数的限制，对超过2维的非线性方程组的求解不太合适，需要寻求其他的求解方法。,一般非线性方程的数值解,MATLAB中有两个求非线性方程数值解的函数：函数fzero()采用数值解法求解非线性方程。zfzero(fun,x0)fsolve()函数则采用非线性最小二乘算法求解非线性方程组。xfsolve(fun,x0)使用这两个函数一般分为两个步骤：第一步定义函数文件，并存为.m文件；第二步调用函数命令求解。,fzero()函数示例,第一步：编写M文件，并存为在文件f.m中。functiony=f(x)y=x.3-2*x-5;第二步：在Command窗口中输入命令解方程。,z=2.0946,z=fzero(f,2),注意：该函数为多项式，实际的零点有3个，其余两个为一对共轭复零点。由于fzero不能求复数零点，因此，任意的初始估计都将给出同样的结果。fzero()函数认为零点是函数穿过x轴的点，对于与x轴相切的点，该函数将不能计算出来。对于非连续函数，fzero()函数返回非连续点的值代替零点。,fsolve()函数示例,验证结果：,第一步：编写M文件，并存为在文件fun.m中。functionF=fun(x)F=2*x(1)-x(2)-exp(-x(1);-x(1)+2*x(2)-exp(-x(2);第二步：在Command窗口中输入命令解方程。,Optimizationterminated:first-orderoptimalityislessthanoptions.TolFun.x=0.56710.5671,x0=-5*ones(2,1);options=optimset(Display,iter);%设置显示输出中间结果x=fsolve(fun,x0,options),f=1.0e-006*-0.4059-0.4059,f=fun(x),迭代中间结果,迭代中间结果,NormofFirst-orderTrust-regionIterationFunc-countf(x)stepoptimalityradius0347071.22.29e+00411612003.415.75e+0031293147.0211.47e+0031312854.45213881415239.5271107151867.0412130.8162116.704219.0517242.4278812.2618270.0326580.7595110.2062.59307.03149e-0060.1119270.002942.510333.29525e-0130.001691326.36e-0072.5Optimizationterminated:first-orderoptimalityislessthanoptions.TolFun.,inline()函数和匿名函数调用法,方法1：inline()函数,方法2：匿名函数,f1=inline(P(1)2+P(2)2-1;0.75*P(1)3-P(2)+0.9,P);x,y=fsolve(f1,1,2),f2=(P)P(1)2+P(2)2-1;0.75*P(1)3-P(2)+0.9;x,y,c,d=fsolve(f2,-1,0);x,y,kk=d.iterations,x=0.35700.9341y=1.0e-009*0.12150.0964,x=-0.98170.1904y=1.0e-010*0.5625-0.4500kk=4,例：绝热连续搅拌釜式反应器的转化率,计算达到稳定状态的温度和转化率？,理论分析,理论分析,理论分析,MATLAB求解方法,作图法（1）,作图法（2）,源程序（1）作图法,%主程式functioncstr_nesclearclc%koCAO=exp(20);ER=10000;tau=0.25;H=250;To=450;%checkforsolution%t=450:700;y1=(t-450)/250;a=.25*exp(20-10000./t);,b=1./a+2;y2=.5*(b-sqrt(b.2-4);plot(t,y1,r,t,y2,c)xlabel(temperature)ylabel(conversionrate)title(massbalancevs.energybalance)%findsolutions:usefsolve(共有三组解),作图法结果,源程序（2）数值求解,%findsolutions:usefsolve(共有三组解)%X1=fsolve(fun_cstr,0.9650);%0.9是指转化率x的初始猜测值;%650是指出口温度T的初始猜测值X2=fsolve(fun_cstr,0.3500);X3=fsolve(fun_cstr,0.1450);%printout%disp(sprintf(nthepossiblesolutionsare:n)disp(conv.temp.)disp(X1)disp(X2)disp(X3),%函数式:fun_cstr.mfunctiony=fun_cstr(X)%传入的X有两个:x,T%globalkoCA