《数学建模电子教案.ppt》由会员分享，可在线阅读，更多相关《数学建模电子教案.ppt（160页珍藏版）》请在三一办公上搜索。
1、数学建模电子教案,重庆邮电大学数理学院沈世云,第一节 MATLAB简介,什么是 MATLAB?MATLAB 能干什么?掌握 MATLAB 应用实例,重庆邮电大学-数学建模,学习 MATLAB,什么是 MATLAB?MATLAB 能干什么?掌握 MATLAB 应用实例,重庆邮电大学-数学建模,什么是 MATLAB?,1.MATLAB 代表MATrix LABoratory它的首创者是美国新墨西哥大学计算机系的系主任Cleve Moler博士，他在教授线性代数课程发现其他语言很不方便，篇构思开发了MATLAB。最初采用FORTRAN语言编写，20世纪80年代后出现了MATLAB的第二版，全部采用C
2、语言编写.1984年Moler博士和一批数学家及软件专家创建了MathWorks公司，专门开发MATLAB。1993年出现了微机版，到2003年是6.5版,重庆邮电大学-数学建模,2.一种演草纸式的科学计算语言3.MATLAB 是一高性能的技术计算语言.强大的数值计算和工程运算功能符号计算功能强大的科学数据可视化能力 多种工具箱,重庆邮电大学-数学建模,MATLAB 能干什么？,MATLAB可以进行：数学计算、算法开发、数据采集建模、仿真、原型 数据分析、开发和可视化科学和工程图形应用程序的开发，包括图形用户界面的创建。MATLAB广泛应用于：数值计算、图形处理、符号运算、数学建模、系统辨识、
3、小波分析、实时控制、动态仿真等领域。,重庆邮电大学-数学建模,掌握 MATLAB,MATLAB的构成：MATLAB开发环境：进行应用研究开发的交互式平台MATLAB 数学与运算函数库：用于科学计算的函数MATLAB 语言：进行应用开发的编程工具图形化开发：二维、三维图形开发的工具应用程序接口(API)：用于与其他预言混编面向专门领域的工具箱：小波工具箱、神经网络工具箱、信号处理工具箱、图像处理工具箱、模糊逻辑工具箱、优化工具箱、鲁棒控制工具箱等几十个不同应用的工具箱。,重庆邮电大学-数学建模,开发环境,包括：命令窗口、图形窗口、编辑窗口、帮助窗口。,重庆邮电大学-数学建模,命令窗口,可在提示符
4、后输入交互式命令 结果会自动的产生例如：,重庆邮电大学-数学建模,图形窗口,在窗口中输入：Plot(1,2,4,9,16,1,2,3,4,5)MATLAB 划出如下图形:,编辑窗口,用来创建和修改M-files(MATLAB 脚本),帮助窗口,重庆邮电大学-数学建模,The MATLAB Language,MATLAB 语言的特点Matlab的基本数据单元是不需指定维数的矩阵。Matlab的所有计算都是通过双精度进行的，在内存中的数都是双精度的。double 是一个双精度浮点数，每个存储的双精度数用64位。char用于存储字符，每个存储的字符用16位。,重庆邮电大学-数学建模,MATLAB的程
5、序构成：,重庆邮电大学-数学建模,常变量及其命名规则,变量名可以有数字、字母、下划线构成；变量的首字符必须是字母；区分变量名的大小写每个变量名最长只能包含19个字符。,重庆邮电大学-数学建模,一.矩阵：1.矩阵的建立与表示法：在命令窗口中输入：A=1,2,3;4,5,6;7,8,9 可以得到：A=1 2 3 4 5 6 7 8 9若要显示整行或整列，则可以用(:)冒号来表示。冒(:)代表矩阵中行(ROWS)或列(COLUMNS)的全部。例如执行命令：A(:,2),就会显示第2列的全部，结果为：ans=2 5 8,重庆邮电大学-数学建模,其他特殊矩阵的生成方法：1）、eye(m,n)或eye(m
6、)产生m*n 或 m*m的单位矩阵。例如：eye(3,4)与eye(3)分别产生如下矩阵：1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 12）、zeros(m,n)或 zeros(m)产生m*n 或m*m 的零矩阵。例如：zeros(3,4)与zeros(3)分别产生如下矩阵：0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,重庆邮电大学-数学建模,3）、ones(m,n)或ones(m)产生m*n或m*m的全部元素为1的矩 阵。例如：ones(3,4)与ones(3)分别产生如下矩阵：1 1 1 1 1 1 1 1 1 1
7、1 1 1 1 1 1 1 1 1 1 12.常用矩阵函数：1）、d=eig(A)返回矩阵A的特征值所组成的列向量；v,d=eig(A)返回分别由矩阵A的特征向量和特征值（以其为主对角线元素，其余元素为零）的两个矩阵。例如执行命令：v,d=eig(A)结果为：,重庆邮电大学-数学建模,v=d=0.2320 0.7858 0.4082 16.1168 0 0 0.5253 0.0868-0.8165 0-1.1168 0 0.8187-0.6123 0.4082 0 0-0.0000 其中v(:,i)为d(i,i)所对应的特征向量。2）、det(A)计算行列式A的值。例如：det(A)结果为：a
8、ns=0,重庆邮电大学-数学建模,3）、expm(A)对矩阵A求幂。例如：expm(A)结果为：ans=1.0e+006*1.1189 1.3748 1.6307 2.5339 3.1134 3.6929 3.9489 4.8520 5.75524）、inv(A)求矩阵A的逆。例如：inv(A)结果为：Warning:Matrix is close to singular or badly scaled.Results may be inaccurate.RCOND=2.055969e-018.ans=1.0e+016*-0.4504 0.9007-0.4504 0.9007-1.8014 0
9、.9007-0.4504 0.9007-0.4504,重庆邮电大学-数学建模,5）、orth(A)返回对应于A的正交化矩阵。例如：orth(A)结果为：ans=0.2148 0.8872 0.5206 0.2496 0.8263-0.38796）、poly(A)若A为一矩阵，则返回A的特征多项式。例如：poly(A)结果为：ans=1.0000-15.0000-18.0000-0.0000 若A为一向量，则返回以A的元素为根的特征多项式。例如：r=1,2,3;p=poly(r)结果为：p=1-6 11-67）、rank(A)计算矩阵A的秩。例如：r=rank(A)结果为：r=2,重庆邮电大学-
10、数学建模,3.矩阵的四则运算符号：加“+”减“”乘“*”除“/”共轭转置“”非共轭转置“.”例如：b=1+2i;3+4ib=1.0000+2.0000i 3.0000+4.0000ibans=1.0000-2.0000i 3.0000-4.0000ib.ans=1.0000+2.0000i 3.0000+4.0000i,重庆邮电大学-数学建模,4.矩阵分解：1）、q,r=qr(A)将矩阵A做正交化分解，使得A=q*r。q为单位矩阵(unitary matrix),其范数（norm）为1。r为对角化的上三角矩阵。例如：q,r=qr(A)q=-0.1231 0.9045 0.4082-0.4924
11、 0.3015-0.8165-0.8616-0.3015 0.4082 r=-8.1240-9.6011-11.0782 0 0.9045 1.8091 0 0-0.0000 norm(q)ans=1.0000,重庆邮电大学-数学建模,2）、L,U=lu(A)将矩阵A做对角线分解，使得A=L*U，L为下三角矩阵（lower triangular matrix），U为上三角矩阵(upper triangular matrix)。例如：L,U=lu(A)L=0.1429 1.0000 0 0.5714 0.5000 1.0000 1.0000 0 0 U=7.0000 8.0000 9.0000
12、0 0.8571 1.7143 0 0 0.0000,重庆邮电大学-数学建模,二、多项式：多项式是用向量形式来表示，从最右边算起，第一个为0阶系数，第二个为1阶系数，依次类推。例如一个一元三次多项式：4x3+3x2+2x+1 用4 3 2 1来表示。1.多项式的运算：1）、乘：conv指令执行多项式的相乘运算，指令格式为：z=conv(x,y)例如：x=1 3 5;y=2 4 6;z=conv(x,y)z=2 10 28 38 30,重庆邮电大学-数学建模,如果要对两个以上的多项式进行相乘，以重复使用conv指令，例如：（x,y同上）conv(conv(x,y),x)ans=2 16 68 1
13、72 284 280 1502）、分解：与1）相反，用deconv指令，其指令格式为：z,r=deconv(x,y)表示x除以 y商为z，余数为r。例如：z,r=deconv(z,x)z=2 4 6 r=0 0 0 0 0,重庆邮电大学-数学建模,3）、求根：roots指令用于求多项式的根。例如：fx=1 3 2;rootoffx=roots(fx)rootoffx=-2-1 4）、polyval(p,x)计算多项式p在x出的值，其中x可以是点或向量或矩阵。例如：p=1-6 11-6；x=1;p1=polyval(p,x)结果为：p1=0 x=1,2,3;p2=polyval(p,x)结果为：
14、p2=0 0 0,重庆邮电大学-数学建模,x=A;p3=polyval(p,x)结果为：p3=0 0 0 6 24 60 120 210 336 5）、polyder(p)求p的微分多项式。例如：p=1-6 11-6;dp=polyder(p)dp=3-12 11,重庆邮电大学-数学建模,6）、r,p,k=residue(x,y)求x/y的部分因式分解。若多项式x,y都没有重根,则可把x/y的比值表示为x/y=r1/(s-p1)+r2/(s-p2)+.+rn/(s-pn)+ks例如 用residue指令求x/(x2+3x+2)的部分因式分解：x=1 0;y=1 3 2;r,p,k=residu
15、e(x,y)r=2-1 p=-2-1 k=,重庆邮电大学-数学建模,当输入三个参数 r,p,k 时，该函数又会生成原来的函数。例如：?x,y=residue(r,p,k)x=1 0 y=1 3 2,重庆邮电大学-数学建模,三 符号变量、符号表达式、抽象函数：函数sym用于生成符号变量和符号表达式，如：x=sym(x)a=sym(alpha)分别创建变量x,alpha f=sym(a*x2+b*x+c)创建变量表达式f，但要注意此式并没有自动创建变量a,b,c,x。可以用函数syms对多个变量同时定义，如：syms a b c x 函数sym也可以用来表示确定的函数，如 f=sym(f(x)生成
16、函数f(x)。,重庆邮电大学-数学建模,四 常见符号计算：1.微分：diff是求微分最常用的函数。其输入参数既可以是函数表达式，也可以是符号矩阵。Diff(f,x,n)表示对f关于x求n阶导数。例如：1).下面程序段将生成表达式sin(ax),并分别对其中的x和a求导。?syms a x?f=sin(a*x);?df=diff(f,x)df=cos(a*x)*a?dfa=diff(f,a)dfa=cos(a*x)*x,重庆邮电大学-数学建模,2）、若输入参数为矩阵，将对矩阵中的每个元素求导。?syms a x?A=-sin(a*x),sin(a*x);cos(a*x),cos(a*x)A=-s
17、in(a*x),sin(a*x)cos(a*x),cos(a*x)?dy=diff(A,x)dy=-cos(a*x)*a,cos(a*x)*a-sin(a*x)*a,-sin(a*x)*a,重庆邮电大学-数学建模,3）、可用函数jacobian来计算Jacobi矩阵。?syms r l f?x=r*cos(l)*cos(f);?y=r*cos(l)*sin(f);?z=r*sin(l);?J=jacobian(x;y;z,r l f)J=cos(l)*cos(f),-r*sin(l)*cos(f),-r*cos(l)*sin(f)cos(l)*sin(f),-r*sin(l)*sin(f),r
18、*cos(l)*cos(f)sin(l),r*cos(l),0,重庆邮电大学-数学建模,2.积分：用函数int来求符号表达式的积分。命令格式为：int(f,r,x0,x1)其中f为所要积分的表达式，r为积分变量，若为定积分，则x0,x1为积分上下限。例：?sym x;?sym k real?f=exp(-(k*x)2)f=exp(-k2*x2)?int(f,x,-inf,inf)ans=signum(k)/k*pi(1/2),重庆邮电大学-数学建模,3.级数求和：函数用于对符号表达式求和。例：?syms k;?s1=symsum(1/k2,1,inf)s1=1/6*pi2,重庆邮电大学-数学建
19、模,4.极限：用函数limit来求表达式的极限。函数limit的常用调用格式：数学表达式 命令格式 Limit(f),或limit(f,x)Limit(f,x,a),或 limit(f,a)Limit(f,x,a,left)Limit(f,x,a,right),重庆邮电大学-数学建模,5.化简：1）、collect(f)将表达式中相同次幂的项合并，也可以再输入一个参数指定以哪个变量的幂次合并。2）、expand(f)将表达式展开。3）、horner(f)将表达式转换为嵌套格式。4）、factor(f)将表达式分解因式，并且分解后的多项式的所有系数都为有理数。5）、simplify(f)利用函数
20、规则对表达式进行化简。,重庆邮电大学-数学建模,重庆邮电大学-数学建模,第二节MATLAB的程序设计,重庆邮电大学-数学建模,1 脚本文件和函数文件,11 M脚本文件,对于一些比较简单的问题，在指令窗中直接输入指令计算。,对于复杂计算，采用脚本文件（Script file）最为合适。,MATLAB只是按文件所写的指令执行。,M脚本文件的特点是：,脚本文件的构成比较简单，只是一串按用户意图排列而成的（包括控制流向指令在内的）MATLAB指令集合。,脚本文件运行后，所产生的所有变量都驻留在 MATLAB基本工作空间（Base workspace）中。只要用户不使用清除指令（clear），MATLA
21、B指令窗不关闭，这些变量将一直保存在基本工作空间中。,M文件有两种形式：脚本文件（Script File）和函数文件（Function File）。这两种文件的扩展名，均为“.m”。,重庆邮电大学-数学建模,1 脚本文件和函数文件（续1）,12 M函数文件,与脚本文件不同，函数文件犹如一个“黑箱”，把一些数据送进并经加工处理，再把结果送出来。,MATLAB提供的函数指令大部分都是由函数文件定义的。,M函数文件的特点是：,从形式上看，与脚本文件不同，函数文件的笫一行总是以“function”引导的“函数申明行”。,从运行上看，与脚本文件运行不同，每当函数文件运行，MATLAB就会专门为它开辟一个
22、临时工作空间，称为函数工作空间（Function workspace）。当执行文件最后一条指令时，就结束该函数文件的运行，同时该临时函数空间及其所有的中间变量就立即被清除。,MATLAB允许使用比“标称数目”较少的输入输出宗量，实现对函数的调用。,重庆邮电大学-数学建模,1 脚本文件和函数文件（续2）,13 M文件的一般结构,由于从结构上看，脚本文件只是比函数文件少一个“函数申明行”，所以只须描述清楚函数文件的结构。,典型 M函数文件的结构如下：,函数申明行：位于函数文件的首行，以关键字 functio 开头，函数名以及函数的输入输出宗量都在这一行被定义。,笫一注释行：紧随函数申明行之后以%开
23、头笫一注释行。该行供lookfor关键词查询和 help在线帮助使用。,在线帮助文本区：笫一注释行及其之后的连续以%开头的所有注释行构成整个在线帮助文本。,编写和修改记录：与在线帮助文本区相隔一个“空”行，也以%开头，标志编写及修改该M文件的作者和日期等。,函数体：为清晰起见，它与前面的注释以“空”行相隔。,重庆邮电大学-数学建模,2 函数调用和参数传递,21 局部变量和全局变量,局部（Local）变量：它存在于函数空间内部的中间变量，产生于该函数的运行过程中，其影响范围也仅限于该函数本身。,全局（Global）变量：通过 global 指令，MATLAB也允许几个不同的函数空间以及基本工作空
24、间共享同一个变量，这种被共享的变量称为全局变量。,22 函数调用,在MATLAB中，调用函数的常用形式是：,输出参数1,输出参数2,=函数名(输入参数1,输入参数2,),函数调用可以嵌套，一个函数可以调用别的函数，甚至调用它自己（递归调用）。,重庆邮电大学-数学建模,2 函数调用和参数传递（续）,23 参数传递,MATLAB在函数调用上有一个与众不同之处：函数所传递的参数具有可调性。,传递参数数目的可调性来源于如下两个MATLAB永久变量：,函数体内的 nargin 给出调用该函数时的输入参数数目。,函数体内的 nargout 给出调用该函数时的输出参数数目。,只要在函数文件中包括这两个变量，
25、就可以知道该函数文件调用时的输入参数和输出参数数目。,值得注意：nargin、nargout 本身都是函数，不是变量，所以用户不能赋值，也不能显示。,“变长度”输入输出宗量：varargin、varrgout。具有接受“任意多输入”、返回“任意多输出”的能力。,跨空间变量传递：evalin。,（参考：circle.m，am1.m）,重庆邮电大学-数学建模,23 MATLAB的程序结构和控制流,231 程序结构,循环结构：MATLAB提供两种循环方式。,顺序结构,分支结构：ifelseend。,forend 循环和while-end循环。,232 程序流控制,常用指令：return，echo，i
26、nput，pause，keyboard，break。,switch-case 结构。,try-catch 结构。,警示指令：error，warning。,重庆邮电大学-数学建模,23 MATLAB的程序结构和控制流（续）,233 图形用户界面（GUI）编程,现代的主流应用程序已经从命令行的交互方式转变为以图形界面为主的交互方式，这主要是由于它给用户带来了操作和控制的方便与灵活性。（面向对象编程）,MATLAB能够以比较简单的方式实现一系列的图形界面功能。通过对控件、菜单属性的设置和 Callback 的编写，就能够满足大多数用户的需求。,控件的 Callback 属性：Callback 属性的
27、取值是字符串，可以是某个M文件名或一小段MATLAB语句。当用户激活控件对象（例如：在控件对象图标上单击鼠标左键）时，应用程序就运行 Callback 属性定义的子程序。,菜单的 Callback 属性：Callback 属性的取值是字符串，可以是某个M文件名或一小段MATLAB语句。当用户激活菜单对象时，若没有子菜单就运行 Callback 属性定义的子程序。若有，先运行 Callback 属性定义的子程序，再显示子菜单。,重庆邮电大学-数学建模,24 M文件的调试,编写 M文件时，错误（Bug）在所难免。错误有两种：语法（Syntax）错误和运行（Run-time）错误。,语法错误是指变量
28、名、函数名的误写，标点符号的缺、漏等。对于这类错误，通常能在运行时发现，终止执行，并给出相应的错误原因以及所在行号。,运行错误是算法本身引起的，发生在运行过程中。相对语法错误而言，运行错误较难处理。尤其是M函数文件，它一旦运行停止，其中间变量被删除一空，错误很难查找。,有两种调试方法：直接调试法和工具调试法。,重庆邮电大学-数学建模,24 M文件的调试（续1）,直接调试法：可以用下面方法发现某些运行错误。,在M文件中，将某些语句后面的分号去掉，迫使M文件输出一些中间计算结果，以便发现可能的错误。,在适当的位置，添加显示某些关键变量值的语句（包括使用 disp 在内）。,利用 echo 指令，使
29、运行时在屏幕上逐行显示文件内容。echo on 能显示M脚本文件；echo FunNsme on 能显示名为FunNsme 的M函数文件。,在原M脚本或函数文件的适当位置，增添指令 keyboard。keyboard 语句可以设置程序的断点。,通过将原M函数文件的函数申明行注释掉，可使一个中间变量难于观察的M函数文件变为一个所有变量都保留在基本工作空间中的M脚本文件。,重庆邮电大学-数学建模,第三节 优化问题解法,重庆邮电大学-数学建模,一、无约束优化问题,重庆邮电大学-数学建模,二、线性规划问题,重庆邮电大学-数学建模,三、非线性规划问题,重庆邮电大学-数学建模,重庆邮电大学-数学建模,重庆
30、邮电大学-数学建模,重庆邮电大学-数学建模,第四节图形处理,重庆邮电大学-数学建模,基 本 绘 图 函 数,1 plot,plot3 建立向量或矩阵的图形,2 Loglog x、y轴都取对数标度建立图形,3 Semilogx x轴用于对数标度，y轴线性标度绘制图形,4 Semilogy y轴用于对数标度，x轴线性标度绘制图形,5 Title 给图形加标题,6 Xlabel，Ylabel 给x,y轴加标记,7 Text 在图形指定的位置上加文本字符串,8 gtext 在鼠标的位置上加文本字符串,9 grid 打开网格线,重庆邮电大学-数学建模,绘 图 入 门,1 plot的基本调用格式,2 对於
31、变化剧烈的函数，可用fplot来进行较精确的绘图,3 对符号函数作图可用ezplot,重庆邮电大学-数学建模,特 殊 二（三）维 绘 图 函 数,1 bar(x,y)(barh(x,y),bar3,bar3h)直方（水平）图,2 comet(x,y)(comet3)建立彗星流动图,3 errorbar(x,y,l,u)图形加上误差范围,4 polar(theta,rho)极坐标图,hist(y,x)向量统计的直方图,其中y为要统计的。当x为标 量时，x指定了统计的区间数；当x为向量时，以该向量中 各元素为中心进行统计，区间数等于x向量的长度。,6 rose(theta)极坐标频数累计柱状图,重
32、庆邮电大学-数学建模,9 fill 实心图,7 stairs(x,y)阶梯图,8 stem(x,y,fill)针状图,11 compass 罗盘图,10 feather 羽毛图,quiver,quiver3 向量场图,通常与contour(),gradient()配合使用.,13 pie,pie3 饼图,重庆邮电大学-数学建模,绘 制 三 维 曲 面,mesh(z)语句给出矩阵Z元素的三消隐图，surf和mesh的用法相似。为了方便测试立体绘图，MATLAB提供了一个peaks函数，可产生一个凹凸有致的曲面，包 含了三个局部极大点及三个局部极小点，要画出此函数的最快方法即是直接键入peaks.
33、三维函数还有meshc(),meshz(),surfc(),surfl(),contourf(),waterfall()等.,重庆邮电大学-数学建模,第五节 解方程,重庆邮电大学-数学建模,一 般 的 代 数 方 程,函数solve用于求解一般代数方程的根，假定S为符号表达式，命令solve(S)求解表达式等于0的根,也可以再输入一个参数指定未知数。例：syms a b c x S=a*x2+b*x+c;solve(S)ans=1/2/a*(-b+(b2-4*a*c)(1/2)1/2/a*(-b-(b2-4*a*c)(1/2)b=solve(S,b)b=-(a*x2+c)/x,重庆邮电大学-数
34、学建模,线 性 方 程 组,线性方程组的求解问题可以表述为：给定两个矩阵A和B，求解满足方程AX=B或XA=B的矩阵X。方程AX=B的解用X=AB或X=inv(A)*B表示；方程XA=B的解用X=B/A或X=B*inv(A)表示。不过斜杠和反斜杠运算符计算更准确，占用内存更小，算得更快。,线 性 微 分 方 程,函数dsolve用于线性常微分方程（组）的符号求解。在方程中用大写字母D表示一次微分,D2，D3分别表示二阶、三阶微分，符号D2y相当于y关于t的二阶导数。函数dsolve 的输出方式 格式 说明y=dsolve(Dyt=y0*y)一个方程，一个输出参数u,v=dsolve(Du=v,
35、Dv=u)两个方程，两个输出 参数S=dsolve(Df=g,Dg=h,Dh=-2*f)方程组的解以S.f S.g S.h 结构数组的形式输出,结 果：u=tg(t-c),解 输入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x),结 果 为:y=3e-2xsin（5x）,解 输入命令：x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z,t)；x=simple(x)%将x化简 y=simple(y)z=simple(z),结 果 为：x=(c1-c2+c3+c2e-3t-c3e-3t)e2
36、t y=-c1e-4t+c2e-4t+c2e-3t-c3e-3t+(c1-c2+c3)e2t z=(-c1e-4t+c2e-4t+c1-c2+c3)e2t,非 线 性 微 分 方 程,t，x=solver（f,ts,x0,options）,1、在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分量形式写成.,2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.,注意:,解:令 y1=x，y2=y1,1、建立m-文件vdp1000.m如下：function dy=vdp1000(t,y)dy=zeros(2,1);dy(1)=y(2);d
37、y(2)=1000*(1-y(1)2)*y(2)-y(1);,2、取t0=0，tf=3000，输入命令：T,Y=ode15s(vdp1000,0 3000,2 0);plot(T,Y(:,1),-),3、结果如图,解 1、建立m-文件rigid.m如下：function dy=rigid(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);,2、取t0=0，tf=12，输入命令：T,Y=ode45(rigid,0 12,0 1 1);plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3
38、),+),3、结果如图,图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线.,第六节回归分析,一：一元线性与非线性回归分析,引例：钢材消费量与国民收入的关系,一元回归模型与回归分析,MATLAB软件实现,简介一元非线性回归模型,实验,为了研究钢材消费量与国民收入之间的关系，在统计年鉴上查得一组历史数据。,引例：钢材消费量与国民收入的关系,试分析预测若1981年到1985年我国国民收入以4.5%的速度递增，钢材消费量将达到什么样的水平？,钢材消费量-试验指标(因变量)Y；国民收入-自变量 x；建立数据拟合函数 y=E（Y|x）=f(x)；作拟合曲线图形分析。,问题分析：,钢材消费
39、量y与国民收入x的散点图,回归分析是研究变量间相关关系的一种统计分析。特点：试验指标（因变量）是随机变量。,图形解释：y=E（Y|x）=f(x),假设：（y=E（Y|x）=f(x)）1）Y是一个正态随机变量，即Y服从正态分 布，并且有方差 D(Y)=2。,2）根据观测值作的散点图，观察出函数f(x)是线性形式还是非线性形式。,回归模型及回归分析,1、一元线性回归模型,或,需要解决的问题：1）在回归模型中如何估计参数a、b和2？,知识介绍,2）模型的假设是否正确？需要检验。,3）利用回归方程对试验指标y进行预测或控制？,参数估计,设观测值为(xi,yi)（i=1,2,n),代入模型中，yi=a+
40、bxi+i,最小二乘法:,回归模型的假设检验,提出问题：,1、相关系数检验,|r|1|r|1，线性相关|r|0，非线性相关,模型：Y=a+bx+,H0的拒绝域为：,2、F-检验法平方和分解公式：,记为,认为线性回归效果好,预测与控制,给定的自变量x0，给出E(y0)的点估计量:,y0的置信度为(1)%的预测区间为：,设y在某个区间(y1,y2)取值时,应如何控制x的取值范围,这样的问题称为控制问题。,小结：,或,模型,1、估计参数a，b，2；2、检验模型正确与否；（即b0）3、预测或控制；,已知数据(xi,yi)(i=1,2,n),如何利用MATLAB软件实现以上的统计计算？,MATLAB软件
41、实现,使用命令regress实现一元线性回归模型的计算,b=regress(Y,X)或 b,bint,r,rint,stats=regress(Y,X,alpha),残差及其置信区间可以用rcoplot(r,rint)画图。,默认值是0.05,引例求解,输入：(hg1.m)x=1097 1284 1502 1394 1303 1555 1917 2051 2111 2286 2311 2003 2435 2625 2948 3155 3372;y=698 872 988 807 738 1025 1316 1539 1561 1765 1762 1960 1902 2013 2446 2736
42、 2825;X=ones(size(x),x,pause c,cint,r,rint,stats=regress(y,X,0.05),pausercoplot(r,rint),输出：c=-460.5282(参数a)0.9840(参数b)cint=-691.8478-229.2085(a的置信区间)0.8779 1.0900(b的置信区间)r=79.1248 69.1244-29.3788-104.1112-83.5709-44.5286-109.7219-18.5724-55.6100-23.8029-51.4019 449.6576-33.4128-109.3651 5.8160 92.13
43、64-32.3827(残差向量)rint=（略）（参见残差分析图）stats=0.9631(R2)391.2713(F)0.0000(P0),预测,x1(1)=3372;(hgy1.m)for i=1:5 x1(i+1)=1.045*x1(i);%未来五年国民收入以4.5%的 速度递增 y1(i+1)=-460.5282+0.9840*x1(i+1);%钢材的预 测值endx1,y1,结果,x1=3372.0 3523.7 3682.3 3848.0 4021.2 4202.1y1=3006.8 3162.9 3325.9 3496.3 3674.4,如果从数据的散点图上发现y与x没有直线关系
44、，又如何计算？,例如，试分析年龄与运动（旋转定向）能力,假设模型,一元多项式回归在matlab 软件中用命令polyfit实现。如前面的例子，具体计算如下：,输入：(phg1.m)x1=17:2:29;x=x1,x1;y=20.48 25.13 26.15 30.0 26.1 20.3 19.35 24.35 28.11 26.3 31.4 26.92 25.7 21.3;p,S=polyfit(x,y,2);p,注意：x，y向量的维数要一致。S是一个数据结构，用于其它函数的计算。,计算y的拟合值：输入：Y,delta=polyconf(p,x,S);Y结果：Y=22.5243 26.0582
45、 27.9896 28.3186 27.0450 24.1689 19.6904 22.5243 26.0582 27.9896 28.3186 27.0450 24.1689 19.6904,拟合效果图：,用polytool(x,y,2)还可以得到一个交互式画面。,在工作空间中，输入yhat，回车，得到预测值。,实验内容,1、确定企业年设备能力与年劳动生产率的关系 某市电子工业公司有14个所属企业，各企业的年设备能力与年劳动生产率统计数据如下表。试分析企业年设备能力与年劳动生产率的关系。若该公司计划新建一个设备能力为9.2千瓦/人的企业，估计劳动生产率将为多少？,一矿脉有13个相邻样本点，人
46、为地设定一个原点，现测得各样本点与原点的距离x，与该样本点处某种金属含量y的一组数据如下：,2、测定某矿脉的金属含量,试建立合适的回归模型。（首先画散点图）,二，多元线性与非线性回归分析,引例：某建材公司的销售量因素分析,多元线性回归模型,MATLAB软件实现,简介多元非线性回归模型,实验,1）了解回归分析的基本原理；2）掌握MATLAB的实现方法；3）练习用回归分析方法解决实际问题；,实验目的,某建材公司对某年20个地区的建材销售量Y(千方)、推销开支、实际帐目数、同类商品竞争数和地区销售潜力分别进行了统计。试分析推销开支、实际帐目数、同类商品竞争数和地区销售潜力对建材销售量的影响作用。试建
47、立回归模型，且分析哪些是主要的影响因素。,引例：某建筑材料公司的销售量因素分析,设：推销开支x1 实际帐目数x2 同类商品竞争数x3 地区销售潜力x4,1111.11111,寻找关系：y=E(Y|x1,x2,x3,x4)=f(x1,x2,x3,x4),模型：,假设：,1、因变量Y是随机变量，并且它服从正态分布；2、f(x1,x2,x3,x4)是线性函数(非线性)；,2、多元线性回归模型,模型要解决的问题可归纳为以下几个方面：1）在回归模型中如何估计参数i(i=0,1,m)和2？2）模型的假设(线性)是否正确？3）判断每个自变量xi(i=1,m)对Y的影响是否显著？4）利用回归方程对试验指标 Y
48、进行预测或控制？,知识介绍,参数估计,求解结果,统计分析,2、,3、残差平方和Q，,由此得2的无偏估计,4、对Y的样本方差S2进行分解,回归模型的假设检验,构造F-统计量及检验H0的拒绝域：,注意：衡量y与x1,x2,xm相关程度的指标可以定义复相关系数R，R的值越接近于1，它们的相关程度越密切。,回归系数的检验,主要判断每个自变量xi对y的影响是否显著。,由此可得,MATLAB软件实现,b=regress(Y,X)或 b,bint,r,rint,stats=regress(Y,X,alpha),1、使用命令regress实现多元线性回归,引例求解：,输入：(jzhui.m)x1=5.5 2.
49、5 8 3 8 6 4 7.5 7;(20维）x2=31 55 67 55 70 40 50 62 59;x3=10 8 12 11 11 9 9;x4=8 6 9 16 8 13 11;y=79.3 200.1 135.8 223.3 195;X=ones(size(x1),x1,x2,x3,x4;b,bint,r,rint,stats=regress(y,X),计算结果：（输出）,b=191.9158-0.7719 3.1725-19.6811-0.4501 0 1 2 3 4bint=103.1071 280.7245（系数的置信区间）r=-6.3045-4.2215 8.4422 23
50、.4625 3.3938rint=（略）stats=0.9034（R2）35.0509（F）0.0000（p）Q=r*r2=Q/(n-2)=537.2092（近似）,残差向量分析图,如何分析四个因素x1,x2,x3,x4 对试验指标Y的作用大小？,使用逐步回归方法。在MATLAB软件中使用以下命令：stepwise(X,y,inmodel,alfha),如上例，输入：X=x1,x2,x3,x4;stepwise(X,y,1,2,3),经过观察，得到各种情况下的均方差对比：,最佳回归方程,范例：某化学反应问题,这是一个非线性回归模型的实例,1、问题 为了研究三种化学元素：氢、n戊烷和异构戊烷与生
