数模中常用的matlab语言和方法.docx

1、数模中常用的matlab语言和方法数模培训一、 矩阵的输入二、 循环语句三、 M文件及M函数四、 优化问题五、 回归分析一 矩阵的输入常用的矩阵输入方法有三种:1.直接输入法 2 下载法 3. 编程输入法下面以例解释下载法例1现要根据瑞士地图计算其国土面积。于是对地图作如下的测量：以西东方向为横轴，以南北方向为纵轴。（选适当的点为原点）将国土最西到最东边界在x轴上的区间划取足够多的分点xi，在每个分点处可测出南北边界点的对应坐标y1 ，y2。用这样的方法得到下表x7.010.513.017.534.040.544.548.056.0y1444547505038303034y2445970729

2、3100110110110x61.068.576.580.591.096.0101.0104.0106.5y1363441454643373328y2117118116118118121124121121x111.5118.0123.5136.5142.0146.0150.0157.0158.0y1326555545250666668y2121122116838182868568根据地图比例知18mm相当于40km，试由上表计算瑞士国土的近似面积。（精确值为41288km2）。解题思路：数据实际上表示了两条曲线，实际上我们要求由两曲线所围成的图形的面积。解此问题的方法是数值积分的方法。具体解时

3、我们遇到两个问题：1。数据如何输入；2。没有现成的命令可用。解：对于第一个问题，我们可把数据考备成M文件（或纯文本文件）。然后，利用数据绘制平面图形。键入load mianji.txta=mianji;A=a(1:3,1:9),a(4:6,1:9),a(7:9,1:9);plot(A(:,1),A(:,2),r,A(:,1),A(:,3),g)接下来可以计算面积。键入：a1=trapz(A(:,1)*40/18,A(:,2)*40/18);a2=trapz(A(:,1)*40/18,A(:,3)*40/18);d=a2-a1d = 4.2414e+004至此，问题可以说得到了解决。之所以说还有

4、问题，是我们觉得误差较大。但计算方法的理论给了我们更精确计算方法。只是MATLAB没有相应的命令。想得到更理想的结果，我们可以自己设计解决问题的方法。（可以编写辛普森数值计算公式的程序，或用拟合的方法求出被积函数，再利用MATLAB的命令quad,quad8）2. 下面以例解释编程输入法若矩阵的各元素之间存有某种关系，那么可利用存有的关系，通过编写程序来建立矩阵；还可以利用MATLAB的矩阵命令来建立矩阵。例：输入矩阵 for i=1:5 % 变量i的步长为1，从1变到5 for j=1:5 % 变量j的步长为1,从1变到5 H(i,j)=1/(i+j-1); % 将值1/(i+j-1)赋予矩

5、阵H的第i行第j列的元素end end循环语句MATLAB也有控制流语句,用于控制程序的流程.主要有for循环、while循环、if和break三种控制语句.虽然语句很少,但功能很强. for循环语句for循环语句的一般表达形式为：for i=表达式 可执行语句1 可执行语句nend例：求S=1+2+3+50,可编程如下：s=0;for k=1:50 s=s+k;end while循环while循环语句用来控制一个或一组语句在某逻辑条件下重复预先确定或不确定的次数.while循环语句的一般表达形式为： while 表达式 循环体语句 end例：对于上面同样的问题,可编程如下： S=0; k=0

6、;while kep k=k+1;x=exp(-x0);dx=abs(x-x0)/(1+abs(x);x0=x;end if和break语句MATLAB中if和break语句的作用与使用方式同其它编程语言一样,用来将控制流程进行分流与中断退出.例：可以把上面的解方程的例子中的循环语句改写成: x0=0;while(1) k=k+1; x=exp(-x0); dx=abs(x-x0)/(1+abs(x); if dx=eps break;end x0=x; end 程序中 while(1)说明循环条件总是真,直到满足dx=eps条件跳出循环体. if else end分支结构分支结构有三种形式：

7、 if 表达式 执行语句 end如果表达式的值非0,则执行下面的语句.否则执行end后面的语句.if 表达式 执行语句1 else 执行语句2 endif 表达式1 执行语句1 elseif 表达式2 执行语句2 elseif 表达式3 执行语句3 else （此句可以省略） 执行语句n end例：对函数可编如下的程序：if x=1&x10y=2*x-1;elsey=3*x-11;end练习:求H矩阵的每一行元素的和,最大值,最小值.求H矩阵的每一行元素的和,最大值,最小值.求H矩阵的主对角线元素的和,最大值,最小值.创建M文件创建M文件是MATLAB中的非常重要的内容.事实上,正是由于在MA

8、TLAB 工具箱中存放着大量的M文件,使得MATLAB在应用起来显得简单、方便,且功能强大. M文件有两种形式：命令文件和函数文件当用户要运行的命令较多时,如果直接在命令窗口中逐条输入和运行,有诸多不便.此时可通过编写命令文件来解决这个问题.另外,从前面的许多例子可以看到：MATLAB的许多命令,需要用户通过编写函数文件来执行. 命令文件的创立进入MATLAB命令窗口后,选择“file”下拉式菜单中的“new”进入编辑/调试器（Editer/Debugger）,在编辑/调试器中,编写符合语法规则的命令.编写完命令文件后,选择“file”下拉式菜单中的“save”项,然后依提示输入一个文件名.至

9、此,完成了命令文件的创建. 函数文件的创立函数文件的创立方法与命令文件的创立方法完全一样,只是函数文件的第一句可执行语句是以function引导的定义语句,并且输入文件名时要与定义语句中的函数名相同.建立了函数文件或命令文件后,只要在命令窗口键入命令文件名或函数名,就可执行M文件中所包含的所有命令.下面分别创建并运行一个命令文件和一个函数文件,以了解M文件的创建和运行的全过程.计算所有小于1000的Fibonnaci数.命令文件的创建和运行： 在MATLAB的命令窗口点击“新建”工具栏或在“file”下拉菜单中选“New”中的“M-file”项,进入编辑/调试器. 在编辑/调试器中,输入以下命

10、令：% 计算小于1000的Fibonnaci数 f=1,1; i=1; while f(i)+f(i+1)1000 f(i+2)=f(i)+f(i+1); i=i+1; end f,i 在“file”下拉菜单中选“Save”项,依提示输入文件名“fibno”至此,完成了命令文件fibno.m的创建. 执行fibno.在MATLAB窗口中输入fibno并敲回车键,计算机依次执行fibno中的各条命令后显示如下的结果：ans = Columns 1 through 12 1 1 2 3 5 8 13 21 34 55 89 144Columns 13 through 16 233 377 610

11、987函数文件的创建和运行： 在MATLAB的命令窗口点击“新建”工具栏或在“file”下拉菜单中选“New”中的“M-file”项,进入编辑/调试器. 在编辑/调试器中,输入以下命令：function f=ffibno(n)f=1,1;i=1;while f(i)+f(i+1)nf(i+2)=f(i)+f(i+1);i=i+1;endf 在“ffile”下拉菜单中选“Save”项,依提示输入文件名“ffibno”至此,完成了函数文件ffibno.m的创建. 执行ffibno.在MATLAB窗口中输入ffibno(1000)并敲回车键即可.将上面的关于H矩阵的运算编成M文件和M函数.6：微分方

12、程数值解（单摆问题）单摆问题的数学模型是在初始角度不大时,问题可以得到很好地解决,但如果初始角较大,此方程无法求出解析解.现问题是当初始角为100和300时,求出其解,画出解的图形进行比较。解:若0较小,则原方程可用来近似.其解析解为(t)= 0cost,. 若不用线性方程来近似,那么有两个模型: 取g=9.8,l=25, 100=0.1745, 300=0.5236.用MATLAB求这两个模型的数值解,先要作如下的处理:令x1=,x2=,则模型变为再编函数文件function xdot=danbai(t,x)xdot=zeros(2,1);xdot(1)=x(2);xdot(2)=-9.8/

13、25*sin(x(1);在命令窗口键入t,x=ode45(danbai,0:0.1:20,0.1745,0);t,y=ode45(danbai,0:0.1:20,0.5236,0);plot(t,x(:,1),r,t,y(:,1),k);参考书:数学实验,高等教育出版社7解优化问题线性规划有约束极小问题模型用命令x, fval= linprog(f,A,b,A1,b1,lb,ub)例1：Find x that minimizes f(x)=-5x1-4x2-6x3subject tox1-x2+x3203x1+2x2+4x3423x1+2x2300x1, 0x2,0x3First, enter

14、 the coefficients:f = -5; -4; -6A = 1 -1 1 3 2 4 3 2 0;b = 20; 42; 30;lb = zeros(3,1);Next, call a linear programming routine:x,fval,exitflag,output,lambda = linprog(f,A,b,lb)；Entering x, fval,lambda.ineqlin, and lambda.lower getsx = 0.0000 15.0000 3.0000fval = -78.0000和其它信息。例2：解问题把问题极小化并将约束标准化键入c=-

15、2,-3,5;a=-2,5,-1;b=-10;a1=1,1,1;b1=7;LB=0,0,0;x，y=linprog（c,a,b,a1,b1,LB）得当X=(6.4286,0.5714,0.0000)时,z=-14.5714最大.例3: 解问题解:键入c=-2,-1,1;a=1,4,-1;2,-2,1;b=4;12;a1=1,1,2;b1=6;lb=0;0;-inf;ub=inf;inf;5;x,z=linprog(c,a,b,a1,b1,lb,ub) 得当X=(4.6667,0.0000,0.6667)时,z=-8.6667最小.非线性规划有约束极小问题模型1模型2：MATLAB求解此问题的命

16、令是：x,fval,exitflag,output,lambda,grad,hessian=fmincon(fun,x0,A,b,A1,b1,LB,UB,nonlcon,options,p1,p2,)fun是目标函数的m_文件名.nonlcon是约束函数C(x)和C1(x)的m_文件名.文件输出为C,C1.例:求解最优化问题第1步: 建立目标函数和非线性约束的m_文件.function y=e1511(x)% 目标函数的m_文件y=exp(x(1)*(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2*x(2)+1);function c1,c2=e1511b(x)% 非线性约束的m_

17、文件c1=1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10;c2=;第2步: 运行程序.键入x0=-1,1;a1=1,1;b1=0;x,f,exitflab,output=fmincon(e1511,x0,a1,b1,e1511b)得结果.输出结果的意义:经过4次迭代(iterations:4)收敛到了(exitfag=1)最优解x(1)=-1.2247,x(2)=1.2247,目标函数最优值为1.8951.8回归分析前面我们曾学过拟合。但从统计的观点看，对拟合问题还需作回归分析。例如：有描述问题甲和问题乙的两组数据(x,y)和(x,z)。设x=0 1 2 3 4 ;

18、y=1.0 1.3 1.5,2.0 2.3;z=0.6,1.95,0.9,2.85,1.8;如果在平面上画出散点图，plot(x,y,r.,x,z,b*)那么问题甲的5个点基本在一条直线上而问题乙的5个点却很散乱。如果都用命令c1=polyfit(x,y,1),c2=polyfit(x,z,1)来拟合，将得到同一条直线。 自然对问题甲的信任程度会高于对问题乙的信任程度。所以有必要对所得结果作科学的评价分析。回归分析就是解决这种问题的科学方法。x=0 1 2 3 4 ;y=1.0 1.3 1.5,2.0 2.3;z=0.6,1.95,0.9,2.85,1.8;x=ones(5,1),x;b,bi

19、nt,r,rint,stats=regress(y,x,0.05);b = 0.9600 0.3300bint = 0.7638 1.1562 0.2499 0.4101输出参数的意义：第一个参数的置信区间：0.7638 1.1562第二个参数的置信区间：0.2499 0.4101stats = 0.9829 171.9474 0.0010 0.0063R2=0.9829 F=171.9474, p=0.0010.R是衡量y与x的相关程度的指标，称为相关系数。R越大，x与y关系越密切。通常R大于0.9才认为相关关系成立。F是一统计指标p是与F对应的概率，当 p0.05时，回归模型成立。此例中

20、p=0 0.00100.05，所以，所得回归模型成立。再输入rcoplot(r,rint)得图形,说明第3个数据应剔除,将会得到更准确的模型.对于第2条曲线做同样的工作.b,bint,r,rint,stats=regress(z,x,0.05);bb = 0.9600 0.3300bint = -1.1142 3.0342 -0.5168 1.1768置信区间大很多.stats = 0.3389 1.5381 0.3031 0.7080说明模型不可靠.下面结合三个具体的例子介绍MATLAB实现回归分析的命令。1 合金强度y与其中含碳量x有密切关系，如下表x0.100.110.120.130.1

21、40.150.160.170.180.200.210.23y42.041.545.045.545.047.549.055.050.055.055.560.5根据此表建立y(x)。并对结果作可信度进行检验、判断x对y影响是否显著、检查数据中有无异常点、由x的取值对y作出预测。解1：在x-y平面上画散点图，直观地知道y与x大致为线性关系。用命令polyfit(x,y,1)可得y=140.6194x+27.0269。作回归分析用命令b,bint,r,rint,ststs=regress(y,x,alpha) 可用help查阅此命令的具体用法残差及置信区间可以用rcoplot(r,rint)画图设回归

22、模型为 y=0+1x,在MATLAB命令窗口中键入下列命令进行回归分析x=0.1:0.01:0.18;x=x,0.2,0.21,0.23;y=42,41.5,45,45.5,45,47.5,49,55,50,55,55.5,60.5;X=ones(12,1),x;b,bint,r,rint,stats=regress(y,X,0.05);b,bint,stats,rcoplot(r,rint)得结果和图b = 27.0269 140.6194bint = 22.3226 31.7313111.7842 169.4546 stats =0.9219 118.0670 0.0000结果含义为0=2

23、7.0269 1=140.61940的置信区间是 22.3226,31.73131的置信区间是 111.7842,169.4546R2=0.9219 F=118.0670, p10-4.R是衡量y与x的相关程度的指标，称为相关系数。R越大，x与y关系越密切。通常R大于0.9才认为相关关系成立。F是一统计指标p是与F对应的概率，当 p0.05时，回归模型成立。此例中 p=0 10-40.05，所以，所得回归模型成立。观察所得残差分布图，看到第8个数据的残差置信区间不含零点，此点视为异常点，剔除后重新计算。此时键入：X(8,:)=;y(8)=;b,bint,r,rint,stats=regress

24、(y,X);b,bint,stats,rcoplot(r,rint)得：b = 27.0992 137.8085bint = 23.8563 30.3421 117.8534 157.7636stats =0.9644 244.0571 0.0000可以看到：置信区间缩小；R2、F变大，所以应采用修改后的结果。2 将17至29岁的运动员每两岁一组分为7组，每组两人测量其旋转定向能力，以考察年龄(x)对这种运动能力(y)的影响。现得到一组数据如下表年龄17192123252729第一人2048251326153002612031935第二人243528112633142692257213试建立关

25、系y(x)，并作必要的统计分析。解2：在x-y平面上画散点图，直观地知道y与x大致为二次函数关系。x=17:2:29;x=x,x;y=20.48,25.13,26.15 30,26.1,20.3,19.35,24.35,28.11,26.3,31.4,26.92,25.7,21.3;plot(x,y,*)设模型为y=a1x2+a2x+a3此问题可以利用命令polyfit(x,y,2)来解，也可以像上题一样求解。polyfit(x,y,2)ans = -0.2003 8.9782 -72.2150即t=17:0.2:29;d=-0.2003*t.2+8.9782*t-72.215;plot(t,

26、d)下面介绍用命令polytool来解。首先在命令窗口键入polytool(x,y,2)得到一个交互式窗口窗口中绿线为拟合曲线、红线为y的置信区间、可通过移动鼠标的十字线或通过在窗口下方输入来设定x值，窗口左边则输出与x对应的y值及y的置信区间。通过左下方的Export下拉菜单可输出回归系数等。更详细的解释可通过help查阅。3 某厂生产的某产品的销售量与竞争对手的价格x1和本厂的价格x2有关。下表是该产品在10个城市的销售记录。x1120140190130155175125145180150x210011090150210150250270300250y(个)1021001207746932

27、6696585试建立关系y(x1，x2)，对结果进行检验。若某城市本厂产品售价160（元），对手售价170（元），预测此产品在该城市的销售量。解3.这是一个多元回归问题。若设回归模型是线性的，即设y=0+1x1+2x2那么依然用regress(y,x,alpha)求回归系数。键入x1=120,140,190,130,155,175,125,145,180,150;x2=100,110,90,150,210,150,250,270,300,250;y=102,100,120,77,46,93,26,69,65,85;x=ones(10,1),x1,x2;b,bint,r,rint,stats=regress(y,x);b,bint,stats, b = 66.5176 0.4139 -0.2698bint = -32.5060 165.5411 -0.2018 1.0296 -0.4611 -0.0785stats = 0.6527 6.5786 0.0247 351.0445 得b = 66.5176 0.4139 -0.2698bint = -32.5060 165.5411 -0.20