环境模拟计算的MATLAB程序设计.ppt

1、MATLAB基础及其在环境科学中的应用,培训内容,环境模拟计算的MATLAB程序设计； 环境统计常见数据分析方法的MATLAB实现； 环境基本模型的MATLAB实现及应用（1）； 环境基本模型的MATLAB实现及应用（2）； 非传统环境数学问题的MATLAB实现； 讨论。,环境模拟计算的MATLAB程序设计,第一讲,一、MATLAB基本数学运算,如何操作MATLAB？ 如何在使用中快速获得帮助？ 其主要数值运算功能有哪些？,一、MATLAB基本数学运算,工作窗口系统 ：经常用到的窗口有6个：命令窗口、M文件编辑窗口、工作空间窗口、当前目录窗口、工具箱窗口、命令记录窗口,一、MATLAB基本数学

2、运算,查询帮助系统：MATLAB的帮助命令有很多，可以分为help系列和lookfor命令两大类。 help系列 包括help、help+函数名、helpwin和helpdesk。 lookfor命令 “help+函数名”的方法对要求已经确切知道函数名，但是由于用户不熟悉的原因，一些函数名可能记忆得不是很清楚，这时候可以用lookfor命令，根据提供的关键字获得相关帮助，具体用法为：lookfor+关键字。需要时间较长。,一、MATLAB基本数学运算,演示帮助系统 进入帮助系统主要有三种方法： 其一就是通过在命令窗输入“demo”命令； 其二是通过Help菜单的“Demos”子菜单； 其三是在

3、命令窗口输入“intro”命令,一、MATLAB基本数学运算,常用操作命令 操作选项（File-优选项）-很有用的！,一、MATLAB基本数学运算,变量与常量 MATLAB变量命名规则： 由字母开头，后面可以跟字母、数字、下划线等；变量名区分大小写； 变量名不能超过31个字符； 许多标点符号在MATLAB中有特殊意义，一般在变量中不允许使用。 变量可以直接赋值，或者依靠计算获得值，然后可以直接调用。,一、MATLAB基本数学运算,变量数据类型 数值型 常用的数值量为双精度浮点数，占8个字节，MATAB表示为doube()。考虑到一些特殊应用，如图像处理，MATLAB还引入了无符号8位整型数据，

4、表示为uint8()，可大大减少其存储空间。此外，MATLAB中还使用其它数据类型，如int8()、int16()、int32() 、int32()。对于MATLAB中数据显示格式，可以用format命令控制，具体用法可以使用help命令查询。 如果需要显示精度较高的数值，可以使用变精度算法函数vpa()来表示。调用格式：vpa(A) 或 vpa(A，n) % A表示数值型变量，n表示显示有效数字位数。 符号型 用于公式推导和一些问题的解析解。符号型变量使用前需要用syms命令申明。语法格式为：syms var_list varform % var_list变量列表；varform变量类型。例

5、如：syms x real y positive %表示定义符号变量x为实型，y为正数。定义多个变量时，变量之间以空格分开，不能用逗号分开。变量类型属可选项，默认为实数变量。 字符串型 MATLAB支持字符串变量，但是要用单引号而不是双引号括起来。,一、MATLAB基本数学运算,语句基本形式 直接赋值语句 基本形式为：变量=赋值表达式 分号“；”的用法 函数文件调用语句 基本形式为：返回变量列表=函数名(输入变量列表) 例如u,v=eig(A)返回矩阵A的特征向量u和特征值v。 命令文件调用语句 调用格式：文件名。 冒号表达式 语句格式：向量=初值：步长：终值限值,v1=0:0.4:pi %注

6、意观察最终值是2.8而不是 v2=0:pi %默认步长为1 v3=0:-0.5:pi %步长为负值，显然不可能生成向量，因此得空矩阵 v4=pi:-0.5:0 %逆序构成向量，步长可以为负值,一、MATLAB基本数学运算,矩阵基本运算 子矩阵提取 语句格式：B=A(S1，S2) %S1、S2分别表示行号、列号构成的向量。 【例】已知矩阵A=1 2 3 4 5;6 7 8 9 0;5 4 3 2 1;0 9 8 7 6，体会下列语句。 A=1 2 3 4 5;6 7 8 9 0;5 4 3 2 1;0 9 8 7 6 %定义矩阵 A1=A(1:2:end, :) %提取矩阵A的1、3行，所有列。

7、关键词end表示最后一行或列 A2=A(3,2,1,1 2 4) %提取矩阵A的3、2、1行，1、2、4列 A3=A(1:2,3:5) %提取矩阵A的1、2行，3、4、5列 A4=A(:,end:-1:1) %提取矩阵A所有行和列，并将矩阵左右翻转，最后一列排到最前 矩阵转置 语句格式：B=A %单引号“”表示对矩阵进行转置操作。 矩阵翻转 语句格式： B=fliplr(A) %左右翻转 B=flipud(A) %上下翻转 B=rot90(A) %逆时针旋转90度,一、MATLAB基本数学运算,矩阵基本运算 矩阵的四则运算 可以直接用“+”、“-”、“*”、“/”、“”表示矩阵之间的四则运算，

8、但要注意矩阵之间运算时维数的相容性。 矩阵乘方 语句格式：B=Ax %A必须是方阵，x可以是正数或分数 点运算 可以直接用“.*”、“./”、“.”分别表示矩阵之间的乘法、除法和乘方的点运算，点运算实际上是对两个参与运算的矩阵中的对应元素的直接运算。,一、MATLAB基本数学运算,关系运算和逻辑运算,一、MATLAB基本数学运算,符号运算:可用于推导数学公式、求解代数方程、积分求解、微分方程的解析解。符号表达式可以通过syms命令、sym函数定义 例如： syms n;h=symsum(3/(2*n-1)2,1,inf) % n从1到无穷大，对表达式3/(2*n-1)2求和 syms x;y1

9、=2*x2-3*x-5;y2=sin(x); fy1=finverse(y1);fy2=finverse(y1); %分别求y1、y2的逆函数 syms x;dsolve(Dx=1+x2,x(0)=2) %D表示微分，微分方程为x=1+x2,初始值小x0=2,一、MATLAB基本数学运算,插值根据由实验和测量得到的由某一函数y=f(x)（可能未知或非常复杂）在不同自变量取值x1、x2xn处的因变量y1、y2yn值，构造一相对简单函数g(x)作为y=f(x)近似表达式：y=f(x)g(x)，使g(x0)=y0，g(x1)=y1g(xn)=yn。g(x)称为插值函数，常采用多项式或分段函数的形式。

10、 一维插值 调用格式：Y=interp1(x,y,X,method) 已知数据(x,y)，用method方法插值，计算确定的插值点X对应得函数值Y。method 插值方法，可取linear（线性插值）、cubic（三次多项式插值）、nearst（临近点插值）、spline（三次样条法插值）。 二维插值 调用格式：Z=interp2(x,y,z,X,Y,method) 已知网格数据(x,y,z)，用method方法插值，计算确定的插值点(X,Y)对应得函数值Z。,X,Y,一、MATLAB基本数学运算,插值 griddata()插值函数 interp2()函数只能处理以网格形式给出的数据。如果已知

11、数据点格式不是以网格形式给出，该函数则无能为力。实际中大部分实测数据不是以网格形式给出的，一般是以多组(xi,yi,zi)的形式给出，为此，MATLAB提供了griddata()函数处理一般的二维插值问题。 调用格式：Z=griddata(x,y,z,X,Y,method) ，式中，x,y,z表示已知的样本点坐标，不要求是网格形式的数据；X，Y表示相应的插值位置(X，Y)，可以是数据点、向量或网格矩阵；method可取linear、cubic、nearest、v4，其中v4是MATLAB4.0版本提供的算法，公认效果较好。,一、MATLAB基本数学运算,已知平面上的数据点可用网格生成函数表示

12、x,y=meshgrid(-5:0.25:5,-5:0.25:5) z坐标可用函数 z=3*(1-x).2.*exp(-(x.2) - (y+1).2) - 10*(x/5 - x.3 - y.5).*exp(-x.2-y.2)表示 现插值求平面网格X,Y=meshgrid(-4:0.2:4,-4:0.2:4)上的点对应的Z 坐标，并绘图。,x,y=meshgrid(-5:0.25:5,-5:0.25:5); z=3*(1-x).2.*exp(-(x.2) - (y+1).2) - 10*(x/5 - x.3 - y.5).*exp(-x.2-y.2); X,Y=meshgrid(-4:0.2

13、:4,-4:0.2:4); Z=interp2(x,y,z,X,Y,cubic); subplot 211; surf(x,y,z); title(原始数据图形); xlabel(X);ylabel(Y);zlabel(Z); subplot 212; surf(X,Y,Z); title(二维插值图形); xlabel(X);ylabel(Y);zlabel(Z);,插值 插值举例,一、MATLAB基本数学运算,插值 例：已知观测数据对(x,y)如下，现求在自变量x0=-1.8 -1.3 -0.7 0.7 1.2 1.7处的y0值。用不同的插值方法，比较其结果。 x= -2.0 -1.5 -

14、1.0 -0.5 0 0.5 1.0 1.5 2.0 y= 0.2 0.31 0.5 0.8 1.0 0.8 0.5 0.31 0.2 【求解】上述问题的MATLAB程序如下。 x=-2.0,-1.5,-1.0,-0.5,0,0.5,1.0,1.5,2.0;%给出自变量观测值 y=0.2,0.31,0.5,0.8,1.0,0.8,0.5,0.31,0.2;%给出因变量观测值 x0=-1.8 -1.3 -0.7 0.7 1.2 1.7; %给出互异插值点 %不同的插值方法 y01=interp1(x,y,x0,linear);y02=interp1(x,y,x0,cubic); y03=inte

15、rp1(x,y,x0,nearst);y04=interp1(x,y,x0,spline); %将插值结果绘图对比 subplot 211;%分割绘图窗为2*1=2个，将绘制图形置于第1个子窗口 plot(x,y,b-,x0,y01,*);%原始观测点画线,插值点画“*” legend(原曲线,插值点);%标注图例 title(linear);%标注标题 xlabel(x轴);%标注x轴 ylabel(y轴);%标注y轴 subplot 212;plot(x,y,b-,x0,y02,*); legend(原曲线,插值点); title(cubic);xlabel(x轴);ylabel(y轴);

16、 cubic和spline的插值结果会更接近理论值，而 linear和nearst的插值结果会比较粗糙。,一、MATLAB基本数学运算,拟合（多项式拟合、线性回归、非线性回归等，此不赘述） 非线性最小二乘拟合 ： k=lsqnonlin(fun,k0) %对函数fun（其值最小前提下），在初始参数为k0情况下，拟合出最佳参数k。,一、MATLAB基本数学运算,拟合-举例 lsqnonlin,已知数据某污染物在反应器中的浓度（c）随反应时间（t）的变化数据如下： t=0.27 0.6 1.0 1.7 2 3 4 6 c=19.27 18.25 16.34 14.30 12.43 9.42 6.4

17、5 5.37。 要求拟合出函数c(t)=c0e-kt中的系数c0和k。,【求解】设c0和k分别用数组K的元素K(1)、K(2)表示，则可以建立M函数，然后调用lsqnonlin函数求解。 首先定义函数c2fun28.m function Zmin=c2fun28(K) global t c; %定义全局变量 t=0.27 0.6 1.0 1.7 2 3 4 6; c=19.27 18.25 16.34 14.30 12.43 9.42 6.45 5.37; Zmin=c-K(1)*exp(-K(2)*t);%定义目标函数:观测值-估计值 然后，调用函数lsqnonlin进行拟合 global

18、t c;%定义全局变量 K0=1,1;%初值均取1 K=lsqnonlin(c2fun28,K0) tt=0:0.5:6;%定义绘图区间 cc=K(1)*exp(-K(2)*tt);%估算拟合值 plot(tt,cc,-rp,t,c,b*); legend(拟合曲线,观测值) 运算结果：K =21.0474 0.2582，c0=21.0474，k=0.2582。,一、MATLAB基本数学运算,规划优化问题,一、MATLAB基本数学运算,规划优化问题-线性规划举例 【求解】给出相关输入参数的值后，调用linprog()函数立即得到结果。 f=-2 -3 -4 -3 -5; A=0 2 1 4 2

19、;3 4 5 -1 -1; B=54;62; Aeq=;Beq=; xl=0 0 3.3 0.7 2.6; xu=; X,fmin,flag,num=linprog(f,A,B,Aeq,Beq,xl,xu,) 运行结果如下： Optimization terminated successfully. X =23.3833 0.0000 3.3000 0.7000 23.9500 fmin =-181.8167;flag = 1 num = iterations: 5 cgiterations: 0 algorithm: lipsol,一、MATLAB基本数学运算,规划优化问题-非线性规划举例,

20、x0=1,1,1;%定义初值 A=;B=;Aeq=8 14 7;Beq=56； xl=0,0,0;xu=;%给出变量值变化范围 X,fmin,flag,num=fmincon (c2fun2_11obj,x0,A,B,Aeq,Beq,xl,xu,c2fun2_11nlin) function Zmin=c2fun2_11obj(x) Zmin=1000-x(1)2-2*x(2)2-x(3)2-x(1)*x(2)-5*x(1)*x(3);%目标函数 function C,Ceq=c2fun2_11nlin(x) C=x(1)*x(2)+x(1)3-30;x(1)*x(3)+x(2)3-40;%非

21、线性不等式约束部分 Ceq=x(1)2+x(2)2+x(3)2-35;%定义非线性等式约束部分 运算结果 Optimization terminated successfully: x =2.1918 0 5.4951；fmin =904.7794；flag =1；num =iterations: 9。,一、MATLAB基本数学运算,概率论与数理统计,正态分布为例：概率密度函数： 概率密度：pn=normpdf(x, , )； 概率分布：F=normcdf(x, , )； 逆概率：x=norminv(F, , ) 其他如泊松分布、二项分布、t分布等不在赘述,例】分别绘制（，）为（-1,1）、（

22、0,0.5）、（0,1）、（0,0.8）、（1,1）时的概率密度函数和概率分布函数曲线。 【求解】首先在-4，4区间构建自变量x（横坐标），再分别定义和2，这样就可以调用相关函数进行计算、绘图。 x=-4:0.2:4;pn=;F=; mu=-1 0 0 0 1;sigma=1 0.5 1 0.8 1; for k=1:length(mu) pn=normpdf(x,mu(k),sigma(k); subplot 211;plot(x,pn);%绘制概率密度函数图 title(概率密度函数); gtext(strcat( ,num2str(mu(k), , ,num2str(sigma(k),

23、); hold on; F=normcdf(x,mu(k),sigma(k); subplot 212;plot(x,F);%绘制概率分布函数图 title(概率分布函数); gtext(strcat( ( , num2str(mu(k), , , num2str(sigma(k), ) ); hold on; end 计算绘图结果如图。从中可以观察：方差相同，则曲线形状相同，只是对不同平均值沿x轴进行平移；若方差不同，则曲线形状不一样，方差越小，曲线越陡，最大概率密度越大。,一、MATLAB基本数学运算,概率论与数理统计问题 统计量分析 均值和中值：m=mean(x)，mm=median(x

24、) 上述函数分别计算向量x的算数均值和中值（中位数）。如果输入参数x为矩阵，则分别计算每列的均值和中值，相当于m=meax(x,1)，mm=median(x,1)。如果要计算每行的均值和中值，则调用格式为m=meax(x,2)，mm=median(x,2)。 极差、标准差、协方差 ：r=range(x)，s=std(x)，c=cov(x) 如果x为向量（行或列均可），则上述函数分别可以计算其极差、标准差、方差。如果输入参数x为矩阵，则分别计算每列的极差、标准差以及列与列之间的协方差（形成协方差矩阵）。 相关系数、峰度、偏度 ：r=corrcoef(x)，s=kurtosis(x)，c=skew

25、ness(x) 对于kurtosis()和skewness()，如果x为矩阵，则其分别计算每列的峰度和偏度，如果x为向量，则只计算出一个峰度或偏度。对于corrcoef()，如果x为矩阵，则其分别计算其列与列之间的相关系数，并形成一个相关系数矩阵；如果输入参数为列向量，则要求要有两个以上的输入列向量，调用格式为r=corrcoef(x, y, z.) 。,二、MATLAB的基本程序设计,M文件 Inline函数 For循环 IF判断 While循环 Switch开关 程序流控制,二、MATLAB的基本程序设计,M文件,M命令文件,若要执行大量的MATLAB命令，可以在M文件编辑窗口中编辑、组织

26、这些命令，然后以扩展名为.m保存在MATLAB目录下的work子目录内（默认为此子目录，也可以保存在其他目录下，运行时用path()函数指定）。 最后在命令窗口调入该文件名即可运行。前面例子中的程序实际上都是这样编制的，因为这样容易调试命令。这就象早期DOS系统中的批处理文件一样。,二、MATLAB的基本程序设计,M文件,M函数文件,利用M函数文件可以将一些经常执行的运算程序编辑为函数的形式，在使用时，直接调入该函数名和输入参数即可，大大减少了重复编写程序代码，精简了程序。M函数文件和M命令文件的显著差别就是M函数文件既可以带输入输出参数，也可以不带，而M命令文件则不能带输入输出参数。,二、M

27、ATLAB的基本程序设计,M文件,M函数文件,二、MATLAB的基本程序设计,应用举例,function RiverCon=c2fun2_30(x) %这是一个计算河流污染物浓度的函数 if nargin1 display(Too many input arguments) end if nargout1 display(Too many output arguments) end M=10000000;u=0.5;A=20;D=50;t=180:60:3600; RiverCon=(M./(A*sqrt(4*D*pi.*t).*exp(-(x-u.*t).2)./(4*D.*t); plot

28、(t,RiverCon, -b*); xlabel(时间); ylabel(下游一定距离处浓度); 这样，只要在命令窗口输入c2fun2_30(500)即可计算500m处浓度过程线。,二、MATLAB的基本程序设计,inline函数,inline()函数的出现实际上是为了克服M函数文件的不足。因为M函数一般需要单独保存，其定义要符合M函数的定义规则等，对于复杂的函数，这样的定义是有必要的。但很多时候函数形式可能比较简单，而且在程序调试的时候，如果被调用函数和主调命令不在一个文件内，会使调试变得繁琐,自变量列表：每个变量用单引号引起来，中间逗号分开,Y = inline(sin(2*pi*X +

29、 sigma), X, sigma),二、MATLAB的基本程序设计,程序设计结构流程,V向量不一定是等间隔的,二、MATLAB的基本程序设计,程序设计结构流程,二、MATLAB的基本程序设计,程序设计结构流程,二、MATLAB的基本程序设计,程序设计结构流程,二、MATLAB的基本程序设计,程序设计结构流程,三、MATLAB的绘图技术,二维绘图 三维绘图 图形编辑,三、MATLAB的绘图技术,基本二维绘图语句,三、MATLAB的绘图技术,基本二维绘图语句,三、MATLAB的绘图技术,基本二维绘图语句,绘制完图形后，还可以用grid on(off)命令在图形上添加（或取消）网格线； 用hold

30、 on(off)命令保持住（或取消保持状态）当前图形，使以后再绘制的图形叠加到原来的图形上（或重新绘图）； 用axis(xmin xmax，ymin ymax)分别指定坐标轴的区限； 用xlabel()、ylabel()命令给坐标轴添加标题； 用title()函数给图形添加标题。,一、二维图形的绘制,其他二维绘图语句,MATLAB还提供了一些特殊意义的绘图函数： loglog(x,y)可以绘制双对数图 semilogx(x,y)可以绘制y轴为对数刻度的二维图形 errorbar(x,y,ymin,ymax)可以在绘制的图形上加上误差范围等,三、MATLAB的绘图技术,二维函数绘图语句,三、MA

31、TLAB的绘图技术,应用举例,三、MATLAB的绘图技术,编程求解,【求解】编制程序如下。 y=inline(200*sin(x)./x,x); subplot(221); fplot(y,1 10);title(Fig1); clear; subplot(222); fplot(tan(x),sin(x),cos(x),2*pi*-1 1 -1 1);title(Fig2);clear; subplot(223); ezplot(x3 + y3 - 5*x*y + 1/5,-3,3,-6 6);title(Fig3); clear syms t ;x= sin(3*t)*cos(t);y=

32、sin(3*t)*sin(t); subplot(224); ezplot(x,y,0,pi); title(Fig4);,一、二维图形的绘制,运行结果,二、三维图形的绘制,基本三维绘图语句,三、MATLAB的绘图技术,基本三维绘图语句,x,y=meshgrid(1:0.5:10,10:0.5:20); z=x.5.*exp(-y.0.5).*y.5.*exp(-x.0.5);,三、MATLAB的绘图技术,基本三维绘图语句,三、MATLAB的绘图技术,三维函数绘图语句,三、MATLAB的绘图技术,三维函数绘图语句,【例】观察下列函数调用后的输出结果。 f =(1-x)2*exp(-(x2)-(

33、y+1)2)- (x/5 - x3 - y5)*exp(-x2-y2)- exp(-(x+1)2 - y2); subplot 221;ezcontour(f,-3,3-6,6,49) subplot 222;ezcontourf(f,-3,3,49) subplot 223;ezsurfc(sqrt(x2 + y2) subplot 224;ff=y/(1 + x2 + y2);ezmeshc(ff,-5,5,-2*pi,2*pi),三、MATLAB的绘图技术,三维函数绘图语句,三、MATLAB的绘图技术,利用绘图函数绘制好图形后，还可以根据需要对图形进行一定的编辑，包括在图形上添加字符、坐标轴名称、标题、箭头，改变线型、颜色，改变图形背景色等等。这种编辑通过三种方式进行： 其一最直接的方式，就是直接点击图形窗口工具栏中相应工具按钮，可以完成相关任务； 其二可以通过图形窗口的“In