MATLAB 求解数学模型的基本知识
目录
1. 熟悉MATLAB 软件运算环境 ................................................................................ 2
（1）运算环境.............................................................................................. 2
（2）怎样定义变量...................................................................................... 2
（3）基本运算.............................................................................................. 3
2. MATLAB中的基本语法 .......................................................................................... 4
（1）矩阵的创建与运算.............................................................................. 4
（2）for 循环语句 . ....................................................................................... 4
（3）if 判定语句 .......................................................................................... 4
3. 函数创建与运行....................................................................................................... 6
4. 数据的读写............................................................................................................... 8
5. 怎样实现函数间数据的传递................................................................................... 9
6. 模型求解（方程组求解）..................................................................................... 10
（1）线性方程组........................................................................................ 10
（2）非线性方程组.................................................................................... 10
7. 实例——传热单元数法预测换热器出口参数..................................................... 12
1. 熟悉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列的零矩阵，矩阵名为X
for i=1:1:10
X(1,i)=i;%给矩阵赋值
end
（3）if 判定语句
例：
syms x y;
x=1;
y=3;
if x==y
disp(' 两者相等！' );
else
disp(' 两者不相等！' );
end
3. 函数创建与运行
MATLAB 与C 语言一样，可以创建函数文件。函数可将复杂程序简化，同时可以重复调用，其基本原理与C 一致。
保存之后，在工作代码区敲入ds(45)，即可将45℃时湿空气的饱和含湿量求解出来。
注意：M 文件的函数名和文件名必须完全一致，否则无法运行。
例一代码：
function f=ds(ta)
%已知空气温度，计算饱和含湿量
ps=2/15*exp(18.5916-3991.11/(ta+233.84));%ps
f=0.622*ps/(101-ps);%ds
end
4. 数据的读写
计算结果需要输出到excel 文档中处理，习惯上采用fopen 函数和fprintf 函数执行数据的输出。
例：
fid=fopen('P001.xls' , 'w' ); %打开文档，w 代表创建‘可写’状态的文档 for i=1:2:100
fprintf(fid,'%.0f\n',i); %把i 输入到P001文档中
end
fclose(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));%ps
f=0.622*ps/(101-ps);%ds
end
function f=ia(ta)
%计算饱和湿空气的比焓值
A=ds(ta);%对ds 函数进行调用，计算出饱和含湿量
f=1.005*ta+A*(2501+1.86*ta);%ia
end
6. 模型求解（方程组求解）
常见的方程组有：线性方程组、非线性方程组、常微分方程组、偏微分方程组等。 常微分方程一般可以转化为非线性方程组，如果不能转化，采用dsolve 函数可以解决一般的常微分方程组。偏微分方程组采用MATLAB 基本上无法求解，需要借助CFD 软件进行求解（对于多元的一阶偏微分可自己离散化进行求解）。 在此只介绍线性方程组、非线性方程组的求解方法。
（1）线性方程组
线性方程可直接采用矩阵计算，MATLAB 求解非常简单。
例如：
⎧x +y +z =1⎪⎨x -y =6
⎪y +z =16⎩
代码：
A=[1 1 1;1 -1 0;0 1 1];%系数矩阵
b=[1;6;16];
X=A\b;
disp(X);%X即为[x,y,z]的矩阵
（2）非线性方程组
很多数学模型都是高维的非线性方程组，计算复杂，变数高。
非线性方程较为复杂，如果方程较少，可采用solve 函数求出所有的解。但当方程较多时，非线性方程有大量的解，solve 函数无法求解。一般采用fsolve 函数进行求解，当然fsolve 函数也可以求解非线性方程。
例：
⎧x +y 2-6=0⎪⎪x ⎨+y -1=0
⎪z
⎪⎩x +y -z =0
代码：
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)对应方程中z
function 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。
控制方程：
NTU =K s A
ξGc p , a
ξ=
ε1=i a '' -i a ' c p , a (t a '' -t a ') 1-exp[-NTU (1-R c )]
1-R c exp[-NTU (1-R c )]
R c =ξGc p , a
M c c p , w
-1 ⎡11⎤K s =⎢+m p n ⎥Av ξBv ⎢c ⎥⎣y 1⎦
G ∆i a =M c ∆t c
…………
代码：
主函数：
function f=danji()
fid=fopen('danji.xls' , 'w' );
x0=[3;80;2;0.5;0.3;33];%赋初值
fprintf(fid,'cse\t Ks\t Rc\t NTU\t E1\t ta1\t 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.3f\t%2.3f\t%2.3f\t%2.3f\t%2.3f\t%2.3f\t%2.3f\n',x,tw1); fprintf('cse\t Ks\t Rc\t NTU\t E1\t ta1\t tw1\n');
fprintf('%2.3f\t %2.3f\t %2.3f\t %2.3f\t %2.3f\t %2.3f %2.3f\n',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);%析湿系数cse
f(2)=1/(1/(21.1*Ua^0.85*x(1)^1.15)+1/(216.6*Vc^0.8))-x(2);%传热系数Ks f(3)=x(1)*G*1.005/(Mc*4.2)-x(3);%热容比Rc
f(4)=x(2)*Ac/(x(1)*G*1.005)-x(4); %传热单元数NTU
f(5)=(1-exp(-x(4)*(1-x(3))))/(1-x(3)*exp(-x(4)*(1-x(3))))-x(5); %热交换效率E1 f(6)=ta0-x(5)*(ta0-TC0)-x(6); %出口气温ta1
end
辅助函数：
function f=ds(TA)
ps=2/15*exp(18.5916-3991.11/(TA+233.84));
f(1)=0.622*ps/(101-ps);
end
function f=ia(TA)
A=ds(TA);
f(1)=1.005*TA+A*(2501+1.86*TA);
end
MATLAB 求解数学模型的基本知识
目录
1. 熟悉MATLAB 软件运算环境 ................................................................................ 2
（1）运算环境.............................................................................................. 2
（2）怎样定义变量...................................................................................... 2
（3）基本运算.............................................................................................. 3
2. MATLAB中的基本语法 .......................................................................................... 4
（1）矩阵的创建与运算.............................................................................. 4
（2）for 循环语句 . ....................................................................................... 4
（3）if 判定语句 .......................................................................................... 4
3. 函数创建与运行....................................................................................................... 6
4. 数据的读写............................................................................................................... 8
5. 怎样实现函数间数据的传递................................................................................... 9
6. 模型求解（方程组求解）..................................................................................... 10
（1）线性方程组........................................................................................ 10
（2）非线性方程组.................................................................................... 10
7. 实例——传热单元数法预测换热器出口参数..................................................... 12
1. 熟悉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列的零矩阵，矩阵名为X
for i=1:1:10
X(1,i)=i;%给矩阵赋值
end
（3）if 判定语句
例：
syms x y;
x=1;
y=3;
if x==y
disp(' 两者相等！' );
else
disp(' 两者不相等！' );
end
3. 函数创建与运行
MATLAB 与C 语言一样，可以创建函数文件。函数可将复杂程序简化，同时可以重复调用，其基本原理与C 一致。
保存之后，在工作代码区敲入ds(45)，即可将45℃时湿空气的饱和含湿量求解出来。
注意：M 文件的函数名和文件名必须完全一致，否则无法运行。
例一代码：
function f=ds(ta)
%已知空气温度，计算饱和含湿量
ps=2/15*exp(18.5916-3991.11/(ta+233.84));%ps
f=0.622*ps/(101-ps);%ds
end
4. 数据的读写
计算结果需要输出到excel 文档中处理，习惯上采用fopen 函数和fprintf 函数执行数据的输出。
例：
fid=fopen('P001.xls' , 'w' ); %打开文档，w 代表创建‘可写’状态的文档 for i=1:2:100
fprintf(fid,'%.0f\n',i); %把i 输入到P001文档中
end
fclose(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));%ps
f=0.622*ps/(101-ps);%ds
end
function f=ia(ta)
%计算饱和湿空气的比焓值
A=ds(ta);%对ds 函数进行调用，计算出饱和含湿量
f=1.005*ta+A*(2501+1.86*ta);%ia
end
6. 模型求解（方程组求解）
常见的方程组有：线性方程组、非线性方程组、常微分方程组、偏微分方程组等。 常微分方程一般可以转化为非线性方程组，如果不能转化，采用dsolve 函数可以解决一般的常微分方程组。偏微分方程组采用MATLAB 基本上无法求解，需要借助CFD 软件进行求解（对于多元的一阶偏微分可自己离散化进行求解）。 在此只介绍线性方程组、非线性方程组的求解方法。
（1）线性方程组
线性方程可直接采用矩阵计算，MATLAB 求解非常简单。
例如：
⎧x +y +z =1⎪⎨x -y =6
⎪y +z =16⎩
代码：
A=[1 1 1;1 -1 0;0 1 1];%系数矩阵
b=[1;6;16];
X=A\b;
disp(X);%X即为[x,y,z]的矩阵
（2）非线性方程组
很多数学模型都是高维的非线性方程组，计算复杂，变数高。
非线性方程较为复杂，如果方程较少，可采用solve 函数求出所有的解。但当方程较多时，非线性方程有大量的解，solve 函数无法求解。一般采用fsolve 函数进行求解，当然fsolve 函数也可以求解非线性方程。
例：
⎧x +y 2-6=0⎪⎪x ⎨+y -1=0
⎪z
⎪⎩x +y -z =0
代码：
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)对应方程中z
function 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。
控制方程：
NTU =K s A
ξGc p , a
ξ=
ε1=i a '' -i a ' c p , a (t a '' -t a ') 1-exp[-NTU (1-R c )]
1-R c exp[-NTU (1-R c )]
R c =ξGc p , a
M c c p , w
-1 ⎡11⎤K s =⎢+m p n ⎥Av ξBv ⎢c ⎥⎣y 1⎦
G ∆i a =M c ∆t c
…………
代码：
主函数：
function f=danji()
fid=fopen('danji.xls' , 'w' );
x0=[3;80;2;0.5;0.3;33];%赋初值
fprintf(fid,'cse\t Ks\t Rc\t NTU\t E1\t ta1\t 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.3f\t%2.3f\t%2.3f\t%2.3f\t%2.3f\t%2.3f\t%2.3f\n',x,tw1); fprintf('cse\t Ks\t Rc\t NTU\t E1\t ta1\t tw1\n');
fprintf('%2.3f\t %2.3f\t %2.3f\t %2.3f\t %2.3f\t %2.3f %2.3f\n',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);%析湿系数cse
f(2)=1/(1/(21.1*Ua^0.85*x(1)^1.15)+1/(216.6*Vc^0.8))-x(2);%传热系数Ks f(3)=x(1)*G*1.005/(Mc*4.2)-x(3);%热容比Rc
f(4)=x(2)*Ac/(x(1)*G*1.005)-x(4); %传热单元数NTU
f(5)=(1-exp(-x(4)*(1-x(3))))/(1-x(3)*exp(-x(4)*(1-x(3))))-x(5); %热交换效率E1 f(6)=ta0-x(5)*(ta0-TC0)-x(6); %出口气温ta1
end
辅助函数：
function f=ds(TA)
ps=2/15*exp(18.5916-3991.11/(TA+233.84));
f(1)=0.622*ps/(101-ps);
end
function f=ia(TA)
A=ds(TA);
f(1)=1.005*TA+A*(2501+1.86*TA);
end