MATLAB求解数学模型的基本知识.docx

MATLAB求解数学模型的基本知识目录1. 熟悉MATLAB软件运算环境2（1）运算环境2（2）怎样定义变量2（3）基本运算32. MATLAB中的基本语法4（1）矩阵的创建与运算4（2）for循环语句4（3）if判定语句43. 函数创建与运行64. 数据的读写85. 怎样实现函数间数据的传递96. 模型求解（方程组求解）10（1）线性方程组10（2）非线性方程组107. 实例传热单元数法预测换热器出口参数121. 熟悉MATLAB软件运算环境（1）运算环境MATLAB帮助文档MATLAB函数的用法在help文档中有详细介绍，看不懂的上网搜集资料。（2）怎样定义变量例：syms x y;x=1;y=exp(x);disp(y);（3）基本运算加减乘除、指数函数、幂函数例：y=x+1;y=x*1;y=1/x;y=exp(x);y=x(-0.5);2. MATLAB中的基本语法MATLAB语法与C语言相近，同时它开发了C语言中没有的大量数学函数，便于科研工作的研究。（1）矩阵的创建与运算例：A=1,2,3;10,20,30;A(2,3);（2）for循环语句例：X=zeros(1,10);%定义一个1行，10列的零矩阵，矩阵名为Xfor i=1:1:10 X(1,i)=i;%给矩阵赋值end（3）if判定语句例：syms x y;x=1;y=3;if x=y disp(两者相等！);else disp(两者不相等！);end3. 函数创建与运行MATLAB与C语言一样，可以创建函数文件。函数可将复杂程序简化，同时可以重复调用，其基本原理与C一致。保存之后，在工作代码区敲入ds(45)，即可将45时湿空气的饱和含湿量求解出来。注意：M文件的函数名和文件名必须完全一致，否则无法运行。例一代码：function f=ds(ta)%已知空气温度，计算饱和含湿量ps=2/15*exp(18.5916-3991.11/(ta+233.84);%psf=0.622*ps/(101-ps);%dsend4. 数据的读写计算结果需要输出到excel文档中处理，习惯上采用fopen函数和fprintf函数执行数据的输出。例：fid=fopen(P001.xls,w);%打开文档，w代表创建可写状态的文档for i=1:2:100 fprintf(fid,%.0fn,i);%把i输入到P001文档中endfclose(fid);%关闭文档若想读取数据，可采用xlsread(filename)读取.xls文档中的数据，或用textread(filename)读取.txt文档中的数据。例：A=xlsread(data001.xls);B=textread(001.txt);A和B均为矩阵。也可以将数据拷贝到MATLAB中，保存为.MAT格式文件，然后再读取。此方法数据稳定性高，一般不会出错，但是读取稍微复杂点。S=load(001.mat);%S是一个结构体，并不是数据矩阵。A=S.data;%结构体S中的data才是矩阵数据（data是.mat数据名）。5. 怎样实现函数间数据的传递有时需要将函数的计算结果，代入到其他函数中进行计算，这时需要交换函数的计算结果。MATLAB传递函数值非常方便，只需直接调用函数名名即可。例：function f=ds(ta)%已知空气温度，计算饱和含湿量ps=2/15*exp(18.5916-3991.11/(ta+233.84);%psf=0.622*ps/(101-ps);%dsendfunction f=ia(ta)%计算饱和湿空气的比焓值A=ds(ta);%对ds函数进行调用，计算出饱和含湿量f=1.005*ta+A*(2501+1.86*ta);%iaend6. 模型求解（方程组求解）常见的方程组有：线性方程组、非线性方程组、常微分方程组、偏微分方程组等。常微分方程一般可以转化为非线性方程组，如果不能转化，采用dsolve函数可以解决一般的常微分方程组。偏微分方程组采用MATLAB基本上无法求解，需要借助CFD软件进行求解（对于多元的一阶偏微分可自己离散化进行求解）。在此只介绍线性方程组、非线性方程组的求解方法。（1）线性方程组线性方程可直接采用矩阵计算，MATLAB求解非常简单。例如：代码：A=1 1 1;1 -1 0;0 1 1;%系数矩阵b=1;6;16;X=Ab;disp(X);%X即为x,y,z的矩阵（2）非线性方程组很多数学模型都是高维的非线性方程组，计算复杂，变数高。非线性方程较为复杂，如果方程较少，可采用solve函数求出所有的解。但当方程较多时，非线性方程有大量的解，solve函数无法求解。一般采用fsolve函数进行求解，当然fsolve函数也可以求解非线性方程。例：代码：function f=xyz (x)%控制方程的函数，将其保存为M函数文件f(1)=x(1)+x(2)2-6;f(2)=x(2)*x(3)+x(2)-1;f(3)=x(1)+x(2)-x(3);end注释：x=x(1),x(2),x(3)，x(1)对应方程中的x，x(2)对应方程中y，x(3)对应方程中zfunction f=m_xyz()%求解方程组的代码，也保存为M函数文件x0=1;1;1;x,fval=fsolve(xyz,x0);disp(x);end 将两个M文件保存后（注意：文件名要与函数相同），在工作代码区输入：m_xyz即可运行函数。7. 实例传热单元数法预测换热器出口参数例一：（单一的换热器，传热单元数法进行迭代求解）已知一个GL换热器的传热系数，进出口参数已知，求其出口参数。进口参数：气温ta0=45，风量200kg/h；冷却水水温TC0=20，冷却水量Mc=600kg/h。控制方程：代码：主函数：function f=danji()fid=fopen(danji.xls,w);x0=3;80;2;0.5;0.3;33;%赋初值fprintf(fid,cset Kst Rct NTUt E1t ta1t tw1n);G=200;%进风质量流量TC0=20;%冷却水初温Mc=420;%冷却水量ta0=45;%入口气温options=optimset(Display,iter);%显示迭代过程x,fval=fsolve(NTU_KES,x0,options,ta0,G,TC0,Mc); %方程求解tw1=G*(ia(ta0)-ia(x(6)/(4.2*Mc)+TC0;%计算出口水温fprintf(fid,%2.3ft%2.3ft%2.3ft%2.3ft%2.3ft%2.3ft%2.3fn,x,tw1);fprintf(cset Kst Rct NTUt E1t ta1t tw1n);fprintf(%2.3ft %2.3ft %2.3ft %2.3ft %2.3ft %2.3f %2.3fn,x,tw1);fclose(fid);end控制方程函数：function f=NTU_KES(x,ta0,G,TC0,Mc)%-参数-Vc=Mc/0.00005/3600/1000;%水流速Ua=G/0.154/3600; %气体流速Ac=14.1; %换热面积%-冷凝器-f(1)=(ia(ta0)-ia(x(6)/1.005/(ta0- x(6)-x(1);%析湿系数csef(2)=1/(1/(21.1*Ua0.85*x(1)1.15)+1/(216.6*Vc0.8)-x(2);%传热系数Ksf(3)=x(1)*G*1.005/(Mc*4.2)-x(3);%热容比Rcf(4)=x(2)*Ac/(x(1)*G*1.005)-x(4); %传热单元数NTUf(5)=(1-exp(-x(4)*(1-x(3)/(1-x(3)*exp(-x(4)*(1-x(3)-x(5); %热交换效率E1f(6)=ta0-x(5)*(ta