数学实验论文

中南财经政法大学期末数学实验设计摘要：数学实验是理工科院校近几十年来的数学教学改革成果。它将经典数学知识，数学建模与计算机应用三者融为一体，使我们能更好地理解数学根本概念，熟悉数学软件。本报告主要是从MATLAB软件操作实验，MATLAB的应用以及MATLAB在数学建模中的应用三个方面来展开。MATLAB软件操作实验主要介绍MATLAB的根本语法和用法，以及它在数学分析〔含微积分〕，高等代数，解析几何，概率论与数理统计中的应用和图形处理功能；MATLAB的应用主要介绍一些简单的实际问题；MATLAB在数学建模中的应用那么是涉及多学科的建模问题，具有一定的使用价值。关键词：MATLAB软件操作实验MATLAB的应用MATLAB在数学建模中的应用一、Matlab简要介绍MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大局部。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言的编辑模式，代表了当今国际科学计算软件的先进水平。MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创立用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。MATLAB的根本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也参加了对C，FORTRAN，C++，JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。20世纪70年代，美国新墨西哥大学计算机科学系主任CleveMoler为了减轻学生编程的负担，用FORTRAN编写了最早的MATLAB。1984年由Little、Moler、SteveBangert合作成立了的MathWorks公司正式把MATLAB推向市场。到20世纪90年代，MATLAB已成为国际控制界的标准计算软件。MATLAB软件的特点：〔1)高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来；〔2)具有完备的图形处理功能,实现计算结果和编程的可视化；〔3)友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握；〔4)功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等),为用户提供了大量方便实用的处理工具。MATLAB优点颇多。MATLAB由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件，其中许多工具采用的是图形用户界面。包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。而且新版本的MATLAB提供了完整的联机查询、帮助系统，极大的方便了用户的使用。简单的编程环境提供了比拟完备的调试系统，程序不必经过编译就可以直接运行，而且能够及时地报告出现的错误及进行出错原因分析。Matlab是一个高级的矩阵/阵列语言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程序〔M文件〕后再一起运行。新版本的MATLAB语言是基于最为流行的C++语言根底上的，因此语法特征与C++语言极为相似，而且更加简单，更加符合科技人员对数学表达式的书写格式。使之更利于非计算机专业的科技人员使用。而且这种语言可移植性好、可拓展性极强，这也是MATLAB能够深入到科学研究及工程计算各个领域的重要原因。MATLAB是一个包含大量计算算法的集合。其拥有600多个工程中要用到的数学运算函数，可以方便的实现用户所需的各种计算功能。函数中所使用的算法都是科研和工程计算中的最新研究成果，而前经过了各种优化和容错处理。在通常情况下，可以用它来代替底层编程语言，如C和C++。在计算要求相同的情况下，使用MATLAB的编程工作量会大大减少。MATLAB的这些函数集包括从最简单最根本的函数到诸如矩阵，特征向量、快速傅立叶变换的复杂函数。函数所能解决的问题其大致包括矩阵运算和线性方程组的求解、微分方程及偏微分方程的组的求解、符号运算、傅立叶变换和数据的统计分析、工程中的优化问题、稀疏矩阵运算、复数的各种运算、三角函数和其他初等数学运算、多维数组操作以及建模动态仿真等。MATLAB自产生之日起就具有方便的数据可视化功能，以将向量和矩阵用图形表现出来，并且可以对图形进行标注和打印。高层次的作图包括二维和三维的可视化、图象处理、动画和表达式作图。可用于科学计算和工程绘图。新版本的MATLAB对整个图形处理功能作了很大的改良和完善，使它不仅在一般数据可视化软件都具有的功能〔例如二维曲线和三维曲面的绘制和处理等〕方面更加完善，而且对于一些其他软件所没有的功能〔例如图形的光照处理、色度处理以及四维数据的表现等〕，MATLAB同样表现了出色的处理能力。同时对一些特殊的可视化要求，例如图形对话等，MATLAB也有相应的功能函数，保证了用户不同层次的要求。另外新版本的MATLAB还着重在图形用户界面〔GUI〕的制作上作了很大的改善，对这方面有特殊要求的用户也可以得到满足。MATLAB对许多专门的领域都开发了功能强大的模块集和工具箱。一般来说，它们都是由特定领域的专家开发的，用户可以直接使用工具箱学习、应用和评估不同的方法而不需要自己编写代码。目前，MATLAB已经把工具箱延伸到了科学研究和工程应用的诸多领域，诸如数据采集、数据库接口、概率统计、样条拟合、优化算法、偏微分方程求解、神经网络、小波分析、信号处理、图像处理、系统辨识、控制系统设计、LMI控制、鲁棒控制、模型预测、模糊逻辑、金融分析、地图工具、非线性控制设计、实时快速原型及半物理仿真、嵌入式系统开发、定点仿真、DSP与通讯、电力系统仿真等，都在工具箱〔Toolbox〕家族中有了自己的一席之地。新版本的MATLAB可以利用MATLAB编译器和C/C++数学库和图形库，将自己的MATLAB程序自动转换为独立于MATLAB运行的C和C++代码。允许用户编写可以和MATLAB进行交互的C或C++语言程序。另外，MATLAB网页效劳程序还容许在Web应用中使用自己的MATLAB数学和图形程序。MATLAB的一个重要特色就是具有一套程序扩展系统和一组称之为工具箱的特殊应用子程序。工具箱是MATLAB函数的子程序库，每一个工具箱都是为某一类学科专业和应用而定制的，主要包括信号处理、控制系统、神经网络、模糊逻辑、小波分析和系统仿真等方面的应用。应用软件开发在开发环境中，使用户更方便地控制多个文件和图形窗口；在编程方面支持了函数嵌套，有条件中断等；在图形化方面，有了更强大的图形标注和处理功能，包括对性对起连接注释等；在输入输出方面，可以直接向Excel和HDF5进行连接。MATLAB系统由MATLAB开发环境、MATLAB数学函数库、MATLAB语言、MATLAB图形处理系统和MATLAB应用程序接口〔API〕五大局部构成。开发环境MATLAB开发环境是一套方便用户使用的MATLAB函数和文件工具集，其中许多工具是图形化用户接口。它是一个集成的用户工作空间，允许用户输入输出数据，并提供了M文件的集成编译和调试环境，包括MATLAB桌面、命令窗口、M文件编辑调试器、MATLAB工作空间和在线帮助文档。数学函数库MATLAB数学函数库包括了大量的计算算法。从根本算法如加法、正弦，到复杂算法如矩阵求逆、快速傅里叶变换等。MATLAB语言是一种高级的基于矩阵/数组的语言，它有程序流控制、函数、数据结构、输入/输出和面向对象编程等特色。图形处理系统图形处理系统使得MATLAB能方便的图形化显示向量和矩阵，而且能对图形添加标注和打印。它包括强大的二维三维图形函数、图像处理和动画显示等函数。应用程序接口MATLAB应用程序接口〔API〕是一个使MATLAB语言能与C、Fortran等其它高级编程语言进行交互的函数库。该函数库的函数通过调用动态链接库〔DLL〕实现与MATLAB文件的数据交换，其主要功能包括在MATLAB中调用C和Fortran程序，以及在MATLAB与其它应用程序间建立客户、效劳器关系。二、Matlab软件根本操作及简单应用〔一〕矩阵的建立和根本运算1.数字变量对于简单的数字运算，可以直接在命令窗口以平常惯用的形式输入。以下是一些关于数字变量的例子：>>a=12,b=345,a*ba=12b=345ans=4140>>r=2,V=4/3*pi*r^3r=2V=33.51032.数组最常用的数组是双精度数值数组。一维数组称为向量，二维数组称为矩阵，一维数组可以视为二维数组的特例。以下是一些关于数组的例子：>>b=[1,3,5,7,9,11]b=1357911>>c=1:2:11c=1357911>>d=linspace(1,11,6)d=1357911>>e=0:5e=012345>>A=[1,2,3;4,5,6;7,8,9]A=123456789>>d=numel(A)d=9>>[n,m]=size(A)n=3m=3>>[i,j]=find(A>3)i=232323j=112233〔二〕符号运算1.字符串与符号型变量在matlab中，所有字符串都用单引号设定后输入或赋值，另外，MATLAB还定义了符号型变量以区别于常规的数值型变量。以下是一些例子：>>x=sym('a')x=a>>y=sym('c*a^2+b-3')y=c*a^2+b–32.符号表达式的假设干运算MATLAB还定义了符号型变量，以区别于常规的数值型变量，可以用于公式推导和数学问题的解析解法，创立符号型变量方法有两种：①用sym(‘’);②用syms。>>symsxy>>S=x^2*y+y^2*x-x^2-x*y+3*yS=x^2*y-x^2+x*y^2-x*y+3*y>>R1=collect(S)R1=(y-1)*x^2+(y^2-y)*x+3*y>>R2=collect(S,y)R2=x*y^2+(x^2-x+3)*y-x^2>>symsxyt>>expand((x-2)*(x-4)*(y-t))ans=8*y-8*t+6*t*x-6*x*y-t*x^2+x^2*y>>symsxy>>F1=factor(x^4-y^4)F1=(x-y)*(x+y)*(x^2+y^2)>>factor(sym('1234567890'))ans=2*3^2*5*3607*3803>>[N,D]=numden(x/y+y/x)N=x^2+y^2D=x*y>>f=1/(1+x^2*y);g=sin(y)g=sin(y)>>C=compose(f,g,x,y)C=1/(y*sin(y)^2+1)>>D=compose(f,g,y,x)D=1/(x^2*sin(x)+1)>>symsx>>f=(1/x^3+6/x^2+12/x+8)^(1/3)f=(12/x+6/x^2+1/x^3+8)^(1/3)>>simple(f)simplify:((2*x+1)^3/x^3)^(1/3)radsimp:(12/x+6/x^2+1/x^3+8)^(1/3)simplify(100):((2*x+1)^3/x^3)^(1/3)combine(sincos):(12/x+6/x^2+1/x^3+8)^(1/3)combine(sinhcosh):(12/x+6/x^2+1/x^3+8)^(1/3)combine(ln):(12/x+6/x^2+1/x^3+8)^(1/3)factor:(12/x+6/x^2+1/x^3+8)^(1/3)expand:(12/x+6/x^2+1/x^3+8)^(1/3)combine:(12/x+6/x^2+1/x^3+8)^(1/3)rewrite(exp):(12/x+6/x^2+1/x^3+8)^(1/3)rewrite(sincos):(12/x+6/x^2+1/x^3+8)^(1/3)rewrite(sinhcosh):(12/x+6/x^2+1/x^3+8)^(1/3)rewrite(tan):(12/x+6/x^2+1/x^3+8)^(1/3)mwcos2sin:(12/x+6/x^2+1/x^3+8)^(1/3)collect(x):(12/x+6/x^2+1/x^3+8)^(1/3)ans=((2*x+1)^3/x^3)^(1/3)>>y=sym('log(x)/sqrt(x)')y=log(x)/x^(1/2)>>dy=diff(y)dy=1/x^(3/2)-log(x)/(2*x^(3/2))>>pretty(dy)1log(x)---------33--22x2x>>symsaxyzt>>S1=findsym(x+i*y-j*z+eps-nan)S1=Emptystring:1-by-0>>S2=findsym(a+t-y,2),S3=findsym(a+t-y,4)S2=y,tS3=y,t,a>>S=a*sin(x)+yS=y+a*sin(x)>>S3=subs(S,x,t),S4=subs(S,x,pi/3)S3=y+a*sin(t)S4=y+(3^(1/2)*a)/23.字符串，符号型变量和数值型变量之间的转换以下是一些字符串，符号型变量和数值型变量之间的转换的例子：>>symsx>>t=1+x;x=1/3x=0.3333>>s=eval(t)s=1.3333>>vpa(s,7)ans=1.333333>>sym('0.3')ans=0.3>>sym(0.3)ans=3/10>>double(s)ans=1.3333〔三〕程序编制及一些常用技巧1.关系运算和逻辑运算MATLAB语言定义了各种比拟关系，关系运算符主要用来比拟数与数、矩阵与矩阵之间的大小，并返回真〔用“1”表示〕、假〔用“0>>2<4ans=1>>4>=7ans=0>>a=4&0a=0>>a=4|0a=1>>a=~4a=0>>a=xor(3,5)a=02.M函数和M文件sin,sum都是MATLAB内嵌的库函数，可以反复调用，十分方便。>>fu=inline('2*x^2+3*x+1')fu=Inlinefunction:fu(x)=2*x^2+3*x+1>>fv=inline('v0+a*t','a','t','v0')fv=Inlinefunction:fv(a,t,v0)=v0+a*t>>v=fv(4,5,2)v=223.MATLAB语言的程序结构与控制语句这里介绍循环结构，转移结构和开关结构。>>fori=1:100s=s+iends=1s=3s=6…s=4950s=5050〔四〕集合和向量的根本运算1.两个集合间的运算>>A=[1,2,3,4;1,2,4,6;6,7,1,4]A=123412466714>>B=[1,2,3,8;1,1,4,6;6,7,1,4]B=123811466714>>C=intersect(A,B,'rows')C=6714>>A=[1,9,6,20];B=[1,2,3,4,6,10,20]B=123461020>>[c,ia,ib]=intersect(A,B)c=1620ia=134ib=157>>A=[1,7,9,6,20];B=[1,2,3,4,6,10,20]B=123461020>>c=setdiff(A,B)c=79>>A=[1,2,3,4]A=1234>>B=[2,4,5,8]B=2458>>C=setxor(A,B)C=1358>>A=[1,2,3,4]A=1234>>B=[2,4,5,8]B=2458>>c=union(A,B)c=1234582.向量的根本运算>>A=[1,2,3]A=123>>B=[3,4,5]B=345>>dot(A,B)ans=26>>a=[1,2,3];b=[4,5,6];c=[7,8,9]c=789>>x=dot(a,cross(b,c))x=0〔五〕多项式和线性方程组的求解1.多项式的表示方式>>p=[1,-4,7,-31]p=1-47-31>>poly2sym(p)ans=x^3-4*x^2+7*x-31>>A=[1,2,3;4,5,6;7,8,0]A=123456780>>p=poly(A)p=1.0000-6.0000-72.0000-27.0000>>poly2sym(p)ans=x^3-6*x^2-72*x-27>>r=[1,3,7]r=137>>p=poly(r)p=1-1131-21>>poly2sym(p)ans=x^3-11*x^2+31*x–212.多项式运算>>p=[2,1,0,-5,3];s=[1,1,1];s0=[0,0,1,1,1]s0=00111>>p+s0ans=211-44>>conv(p,s)ans=233-4-2-23>>[q,r]=deconv(p,s)q=2-1-1r=000-34>>p=[1,3,7]p=137>>roots(p)ans=-1.5000+2.1794i-1.5000-2.1794i>>p=[1,0,3,5,1]p=10351>>d=polyder(p)d=4065>>poly2sym(d)ans=4*x^3+6*x+5>>x=[1,5,7]x=157>>polyval(p,x)ans=1072625843.线性方程组的求解>>A=[1234;4321;1324;4132]A=1234432113244132>>b=[5;4;3;2]b=5432>>x=inv(A)*bx=-1.80001.86673.8667-2.1333>>x1=inv(sym(A))*bx1=-9/528/1558/15-32/15>>norm(double(A*x1-b))ans=0〔六〕图形的绘制1.二维图形绘制例：试绘制出函数y=sin(tanx)-tan(sinx)与y=x+1在x∈[-∏,∏]区间内的曲线。>>x=-pi:0.05:pi>>y=sin(tan(x))-tan(sin(x))>>plot(x,y)>>holdon>>plot(x,x+1,'--r')>>holdoff例：将正弦曲线0~∏/2局部与x轴包围的封闭图形填充成蓝色。>>x=0:pi/60:2*pi;y=sin(x)>>x1=0:pi/60:pi/2;y1=sin(x1)>>plot(x,y,'-r')>>holdon>>fill([x1,pi/2],[y1,0],'b')2.三维图形绘制>>t=0:pi/60:10*pi>>x=sin(t);y=cos(t)>>plot3(x,y,t,'*b')例：绘制带有底座的马鞍面：z=>>x=-8:8>>y=-8:8>>[X,Y]=meshgrid(x,y)>>Z=(X.^2/4^2-Y.^2/5^2)>>meshz(X,Y,Z)例：把三维柱坐标转换成直角坐标。>>theta=0:pi/30:2*pi;ro=sin(theta);[t,r]=meshgrid(theta,ro);>>z=r.*t>>[x,y,z]=pol2cart(t,r,z)>>mesh(x,y,z)〔七〕微积分根本运算1.函数的极限例：求解极限：>>symsxab>>f=x*(1+a/x)^x*sin(b/x)f=x*sin(b/x)*(a/x+1)^x>>limit(f,x,inf)ans=b*exp(a)2.函数的导数例：给定函数,试求出。>>symsxy>>f=y^2*sin(x^2)f=y^2*sin(x^2)>>f1=diff(f,x,2)f1=2*y^2*cos(x^2)-4*x^2*y^2*sin(x^2)>>f2=diff(diff(f,x),y,2)f2=4*x*cos(x^2)3.函数的积分例：求解不定积分与无穷积分。>>symsx>>f1=2*x/(x^3+1)^2;f2=sin(x)/xf2=sin(x)/x>>int(f1,x),int(f2,x,0,inf)ans=(2*x^2)/(3*(x^3+1))-log(x-1/2-(3^(1/2)*i)/2)*(-1/9+(3^(1/2)*i)/9)+log(x-1/2+(3^(1/2)*i)/2)*(1/9+(3^(1/2)*i)/9)-(2*log(x+1))/9ans=pi/24.函数的级数展开例：设试求出该函数在x=0与x=1处的泰勒级数展开的前7项。>>symsx;f=1/(x-2)f=1/(x-2)>>taylor(f,x,7),taylor(f,x,7,1)ans=-x^6/128-x^5/64-x^4/32-x^3/16-x^2/8-x/4-1/2ans=-x-(x-1)^2-(x-1)^3-(x-1)^4-(x-1)^5-(x-1)^6〔八〕代数方程与常微分方程的求解1.代数方程的求解例：用图解法求解方程：。>>ezplot('exp(x)-x^2-10',[0,5])>>holdon,line([0,5],[0,0])2.一般非线性方程数值解例：先用图解法观察方程在区间[0,10]内有多少解，然后试用数值方法求之。>>ezplot('exp(x)-x^2-10',[0,5])>>holdon,line([0,5],[0,0])>>ezplot('5*x^2*sin(x)-exp(-x)',[0,10])>>holdon>>line([0,10],[0,0])>>fun=inline('5*x.^2.*sin(x)-exp(-x)')fun=Inlinefunction:fun(x)=5*x.^2.*sin(x)-exp(-x)>>fsolve(fun,[0,3,6,9])ans=0.50183.14076.28329.42483.常微分方程的符号解例：求解微分方程的通解及在初始条件下的特解。>>dsolve('D2y-Dy-exp(x)=0','x')ans=C3*exp(x)+exp(x)*(x+C2/exp(x))>>dsolve('D2y-Dy-exp(x)-0','y(0)=1','Dy(0)=2','x')ans=exp(x)+x*exp(x)4.常微分方程的数值解法例：求解微分方程。>>fun=inline('-2*y+2*x^2+2*x','x','y')>>[x,y]=ode45(fun,[0,0.5],1)>>plot(x,y)〔九〕优化问题1.无约束最优化问题的数值解法例：二元函数，试用MATLAB提供的求解函数求出其最小值。>>f=inline('(x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2))','x')>>x0=[0;0];ff=optimset;ff.Display='iter'>>x=fminsearch(f,x0,ff)IterationFunc-countminf(x)Procedure01013-0.000499937initialsimplex24-0.000499937reflect36-0.00149944expand47-0.00149944reflect59-0.00349693expand……64121-0.641424contractinside65123-0.641424contractinside66125-0.641424reflect67127-0.641424contractinside68129-0.641424contractinside69131-0.641424contractinside70133-0.641424contractinside71135-0.641424contractinside72137-0.641424contractoutsidex=0.6111-0.30562.有约束最优化问题的求解例：试求解下面的线性规划问题。>>f=[-2,-1,4,-1,1];A=[02142;3-41-12];B=[20;31]B=2031>>Ae=[];Be=[];xm=[0,0,3,1,2]';xM=[]xM=[]>>ff=optimset;ff.LargeScale='off';ff.TolX=1e-15;ff.TolFun=1e-20;ff.Display='iter'>>[x,f_opt,key,c]=linprog(f,A,B,Ae,Be,xm,xM,[],ff)Optimizationterminated.x=14.33334.50003.00001.00002.0000f_opt=-20.1667key=1c=iterations:5constrviolation:8.8818e-016algorithm:'medium-scale:active-set'cgiterations:[]message:'Optimizationterminated.'firstorderopt:4.7370e-015三、数据插值与曲线拟合运算以及数据的根本统计分析〔一〕数据插值与曲线拟合运算1.数据插值〔1〕一维数据插值问题一组样本点(x1,y1),(x2,y2),(x3,y3),…(xn,yn),假设这些样本点全部在某一函数曲线y=f(x)上，那么由这些样本点的信息获得该函数在其他点上的函数值的方法称为函数的插值。以下是一个运用MATLAB语言中interp1函数解决插值问题的例子。year=1900:10:2010;product=[75.995,91.975,105.711,123.203,131.669,150.697,179.323,203.212,226.505,249.633,256.344,267.893];p2005=interp1(year,product,2005)x=1900:1:2010;y=interp1(year,product,x,'spline');plot(year,product,'o',x,y)结果为：p2005=262.1185〔2〕二维数据插值问题以下是一个运用MATLAB语言中interp2函数解决二维插值问题的例子。x=1:6;y=1:4;t=[12,10,11,11,13,15;16,22,28,35,27,20;18,21,26,32,28,25;20,25,30,33,32,30];subplot(1,2,1);mesh(x,y,t)x1=1:0.1:6;y1=1:0.1:4;[x2,y2]=meshgrid(x1,y1);t1=interp2(x,y,t,x2,y2,'spline');subplot(1,2,2);mesh(x1,y1,t1)结果为：2.曲线拟合〔1〕多项式拟合一般多项式拟合的目标是找出一组多项式系数使得多项式能够较好的拟合原始数据。以下是一个运用MATLAB语言中pilyval函数进行拟合的例子。x=[0.5,1.0,1.5,2.0,2.5,3.0];y=[1.75,2.45,3.81,4.80,7.00,8.60];p=polyfit(x,y,3)x1=0.5:0.05:3.0;y1=polyval(p,x1);plot(x,y,'*r',x1,y1,'-b')结果为：p=-0.11561.1681-0.08711.5200〔2〕最小二乘曲线拟合以下是一个运用MATLAB语言中lsqcurvefit函数进行最小二乘拟合的例子。x=0.1:0.1:1;y=[2.3201,2.6470,2.9707,3.2885,3.6008,3.9090,4.2147,4.5191,4.8232,5.1275];Fun=inline('a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4)','a','x');[a,J]=lsqcurvefit(Fun,[1,2,2,3],x,y)结果为：a=2.45872.44891.44662.0734J=0.0018〔二〕数据的根本统计分析1.概率分布及概率问题的求解以下是一些运用MATLAB语言中函数进行进行概率分布及概率问题求解的例子。>>normcdf(1,2,0.5)-normcdf(0,2,0.5)ans=0.0227>>x=norminv(0.6827,2,0.5)x=2.2376>>p=normspec([-2.5,3.5],2,0.5)p=0.99872.统计量分析〔1〕随机变量的均值与方差以下是一个运用MATLAB语言中Rayleigh函数求均值和方差的例子。>>[m,s]=raylstat(0.45)m=0.5640s=0.0869〔2〕随机变量的矩，峰度系数，偏度系数，变异系数以下是一些运用MATLAB语言中函数求随机变量的矩，峰度系数，偏度系数，变异系数的例子。>>x=[9.6,18.5,49.4,23.9,8.1,14,40.7,31.2,33.6,36.8,11.8,9,35,9.4,15];y=[42,54,60,42,1.2,39,68,33.4,53.3,40,31,5.5,46,5.3,23];skewness(x)ans=0.4624kurtosis(x)ans=1.8547skewness(y)ans=-0.4215kurtosis(y)ans=2.25063.大样本数据的处理通常大样本数据处理有以下的处理方法：〔1〕标准化将某矩阵的各行〔列〕元素减去该行〔列〕的均值，再比上该列的标准差。以下是一个运用MATLAB语言中函数将大样本数据标准化的例子。〔2〕都与某一数据进行比拟〔比方都与全国的各项指标进行比拟〕4.直方图与概率值检验函数为了直观地了解随机变量的分布特征，如对称性，峰值等，直方图是广泛使用的方法。以下是一个运用MATLAB语言中函数求直方图的例子。>>a=[459,362,624,542,509,584,433,748,815,505,612,452,434,982,640,742,565,706,593,680,926,653,164,487,734,608,428,1153,593,844,527,552,513,781,474,388,824,538,862,659,775,859,755,649,697,515,628,954,771,609,402,960,885,610,292,837,473,677,358,638,699,634,555,570,84,416,606,1062,484,120,447,654,564,339,280,246,687,539,790,581,621,724,531,512,577,496,468,499,544,645,764,558,378,765,666,763,217,715,310,851]>>hist(a(,30)>>normplot(a()>>mean(a),std(a)ans=600ans=196.6292四、MATLAB在实际问题中的应用〔一〕炮弹发射的角度问题将炮弹发射视为斜抛运动，初始速度为200米/秒，问：要击中水平距离360米，垂直距离160米的目标，当忽略空气阻力时发射角应为多少？进一步思考：如果要考虑水平方向的阻力且设阻力与〔水平方向〕速度成正比，系数是0.1，结果又如何？首先忽略空气阻力建立方程：与。其中，参数有如下方法可以解决：1.图解法利用消元法，消去参变量t，得到。利用matlab软件得到如下结果：>>a=0:0.1:2*pi;y1=160y1=160>>t=1.8./cos(a);y=360*tan(a)-15.876./cos((a).^2)y=1.0e+004*>>plot(a,y,a,y1,’r’,a,t,’b’),grid将图形放大后可得到答案：2.解析法将方程变形得：，由所列方程直接通过matlab迭代可得:>>x=0x=0>>fork=1:20x=atan(0.444+0.0441/(cos(x)^2))endx=0.4629如果要考虑水平方向的阻力且设阻力与水平方向速度成正比，系数为0.1，那么与第一个模型相比主要区别在于水平方向上，炮弹变为变速运动。水平方向上，炮弹运动表示为,即为，初始条件为。对方程进行一系列变换得到：利用matlab来运算得到如下结果：M文件为：functiony=diedai(x)y=160+2000*tan(x)*log(1-36*cos(x)/200)+490/(cos(x))^2*(log(1-36*cos(x)/200))^2>>x0=0;[x,fval]=fzero(@diedai,x0)x=0.473fval=-9.2371e-014>>t=-10/cos(x)*log(1-36*cos(x)/200)t=1.9617〔二〕国土面积的数值计算欧洲一个国家的地图，为了计算出它的国土面积，首先对地图作如下测量：以由西向东方向为x轴，由南向北方向为y轴，选择方便的原点，并将从最西边界点到最东边界点在x轴上的区间适当地分为假设干段，在每个分点的y方向测出南边界点和北边界点的y坐标y1和y2，这样就得到表2.2中的数据〔略〕。采用三次多项式插值计算程序如下：>>x=[7.010.513.017.534.040.544.548.056.061.068.576.580.591.096.0101.0104.0106.5111.5118.0123.5136.5142.0146.0150.0157.0158.0]>>y1=[444547505038303034363441454643373328326555545250666668]>>y2=[4459707293100110110110117118116118118121124121121121122116838182868568]>>newx=7:0.1:158>>L=length(newx)>>newy1=interp1(x,y1,newx,’linear’);newy2=interp1(x,y2,newx,’linear’);>>Area=sum(newy2-newy1)*0.1/18^2*1600Area=4.2414e+004>>fill([newxnewx(L-1:-1:2)],[newy1newy2(L-1:-1:2)],’green’)>>holdon>>plot([x,x],[y1,y2],’*r’)>>xlabel(‘东西距离〔单位：mm〕’),ylabel(‘南北距离〔单位：mm〕’)>>title(‘国土面积计算—三次插值〔比例尺为9:20000000〕’)〔三〕圆周率的计算问题1.迭代加速算法〔1〕程序1functiony=pai1(epsi)%NumericalMethodofPibyusingk=0;n=1;l=1;xn=1;s=1;error=1;epsi=0.0001;whileerror>epsil=-l;n=n+2;xn=1/n;s=s+l*xn;error=abs(4*s-pi);k=k+1;ends=4*sk结果为s=3.1415k=9999〔2〕程序2functiony=pai2(epsi)k=2;n=5;l=1;error=1;epsi=0.00001;s=1-1/3+1/5;whileerror>epsiz=s-1*(n-2)/(4*k*n);l=-l;n=n+2;s=s+l/n;error=abs(4*z-pi);k=k+1;ends=4*zk结果为s=3.1416k=312.公式法functiony=pai3(epsi)k1=0;n=1;s1=1/2;f=1;epsi=0.0001;an=1/2;sn=an;whilesn>epsif=-f;n=n+2;an=an/4;sn=an/n;s1=s1+f*sn;k1=k1+1;endk2=0;n=1;s2=1/3;f=1;an=1/3;sn=an;whilesn>epsif=-f;n=n+2;an=an/9;sn=an/n;s2=s2+f*sn;k2=k2+1;ends=4*(s1+s2)k=k1+k2结果为s=3.1415k=83.几何算法functiony=pai4(epsi)k=0;error=0.5;epsi=0.0000001;c=1/2;L=3*sqrt(3)/2;whileerror>epsic=sqrt((1+c)/2);L=L/c;error=pi-L;k=k+1;endLk结果为L=3.1416k=124.概率的算法functionpi=pai5(N)N=1000;x=rand(N,1);y=rand(N,1);m=sum(x.^2+y.^2<=1);pi=4*m/N;结果为ans=3.0840五、MATLAB在数学建模综合实验中的应用〔一〕施肥效果分析某地区作物生长所需的营养素主要是氮，钾，磷。某作物研究所在该地区对土豆与生菜做了一定数量的实验，实验数据略。1.问题的分析与模型的建立〔1〕做出散点图〔2〕建立回归曲线方程设y与n的函数关系为：。利用最小二乘法我们求下式成立的函数y：。利用MATLAB软件，编程如下：>>n=[03467101135202259336404471]>>y=[15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75]>>polyfit(n,y,2)ans=-0.00030.197114.7416我们得到产量y与n的函数为：同理得到y与p的函数关系为：利用MATLAB软件，编程如下：>>p=[024497398147196245294342]>>y=[33.4632.4736.0637.9641.0440.0941.2642.1740.3642.73]>>polyfit(p,y,2)ans=-0.00010.071932.9161从钾肥的使用量k与产量的关系来看，符合指数函数。其中，a,b,c为待定系数。利用MATLAB软件，编程如下：>>k=[04793140186279372465558651]>>y=[18.9827.3534.8638.5238.4437.7338.4343.8742.7746.22]>>b0=[400.50.05]>>b1=nlinfit(k,y,@examplef,b0)b1=42.66400.54830.0099函数examplef的M-file是functiony=examplef(beta0,x)a=beta0(1);b=beta0(2);c=beta0(3);y=a*(1-b*exp(-c*x));即土豆产量与钾肥之间的函数关系为：。2．模型的应用由于当一种肥料的施肥量变化时，另外的两种肥料的施肥量保持在第7个水平上。于是我们有3种方案：〔n,196，,32〕，〔259，平，372〕，〔259,196，k〕。利用MATLAB软件，编程如下：>>y=inline('-34160*(1-0.56*exp(-0.01*x))+0.64*x+153.37')y=Inlinefunction:y(x)=-34160*(1-0.56*exp(-0.01*x))+0.64*x+153.37>>[x,fval]=fminbnd(y,0,651)x=570.0109fval=-3.3578e+004综上所述，最正确施肥方案为第一个方案〔327.77,196.372〕。3．模型的改良根据农作物的生长原理，氮磷钾3种肥料缺一不可，但同时又是一个有机的整体。因此，要得到农作物的产量与3种肥料之间的使用量的关系，必须考虑3种肥料间的交互影响的数据，也就是说在设计实验时应当采取正交实验，或均匀设计的方法，利用这样实验得到的数据建立农作物的产量与3种肥料之间的多元函数关系，才能比拟准确地找到最正确施肥量。〔二〕最优化捕鱼策略问题试建立数学模型分析如何实现可持续捕捞〔即每年开始捕捞时渔场中各年龄组鱼群条数不变〕，并且在此前提下得到最高的年收获量。1.问题分析表示在t时刻i龄鱼条数k表示4龄鱼捕捞强度n表示每年产卵量表示每年初i龄鱼的数量2．MATLAB模拟程序与结果先建立两个M文件：[buyu.m]functiony=buyu(x)globala10a20asymska10;x1=dsolve('Dx1=-0.8*x1','x1(0)=a10');t=1;a20=subs(x1);x2=dsolve('Dx2=-0.8*x2','x2(0)=a20');t=1;a30=subs(x2);x31=dsolve('Dx31=-(0.8+0.42*k)*x31','x31(0)=a30');t=2/3;a31=subs(x31);x32=dsolve('Dx32=-0.8*x32','x32(2/3)=a31');t=1;a40=subs(x32);x41=dsolve('Dx41=-(0.8+k)*x41','x41(0)=a40');t=2/3;a41=subs(x41);x42=dsolve('Dx42=-(0.8+k)*x42','x42(2/3)=a41');nn=1.109*10^5*(0.5*a31+a41);eq1=a10-nn*1.22*10^11/(1.22*10^11+nn);S=solve(eq1,a10);a10=S(2);symst;t3=subs(subs(int(0.42*k*x31,t,0,2/3)));t4=subs(subs(int(k*x41,t,0,2/3)));total=17.86*t3+22.99*t4;k=x;y=subs(-total);[buyu1.m]globala10a20a[k,mtotal]=fminbnd('buyu',16,18);ezplot(total,0,25)xlabel('捕捞强度系数')ylabel('总收获量〔克〕')title('捕捞强度-总收获量曲线图')formatlong;ktotal=-mtotala10=eval(a10)a20=eval(a20)a30=eval(a30)a40=eval(a40)formatshortclear在MATLAB中输入命令：buyu1，回车得：>>buyu1k=17.362932626767943total=3.887075517793264e+011a10=a20=5.373946224506130e+010a30=a40=8.395506708076184e+007当k=17.36293262688062时，最高年收获为total=3.887075517793262e+011〔克〕此时每年年初四种龄鱼的数量分别为：1.2.5.373946224506130e+0103.4.8.395506708076184e+007〔三〕飞行管理问题在约10000米高空的某边长160公里的正方形区域内，经常有假设干架飞机作水平飞行。区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理。当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区域内的其他飞机发生相撞。如果发生相撞，那么应计算如何调整各架〔包括新进入的〕飞机的飞行方向角，以防止碰撞。1.符号说明t表示时间xi表示第i架飞机的横坐标yi表示第i架飞机的纵坐标dij(t)表示t时刻第i架飞机与第j架飞机间的距离v表示飞机的飞行高度〔v=800〕2.模型的建立由题意可知，目标函数为：约束条件为利用MATLAB来求解，首先写两个M函数：M函数airfun.m：functionf=airfun(delta)f=delta*delta’;M函数airfunco.m：function[c,ceq]=airfunco(delta)x0=[150851501451300];y0=[140,85,155,50,150,0];alpha0=[243236220.515923052]*pi/180;v=800;co=cos(alpha0+delta);si=sin(alpha0+delta);fori=2:6forj=1:i-1t(i,j)=(x0(i)-x0(j))*(co(i)-co(j));t(i,j)=t(i,j)+(y0(i)-y0(j))*(si(i)-si(j));t(i,j)=-t(i,j)/v;t(i,j)=t(i,j)/(co(i)-co(j))^2+(si(i)-si(j))^2;ift(i,j)<0,d(i,j)=1000;elsed(i,j)=(x0(i)-x0(j)+v*t(i,j)*(co(i)-co(j)))^2;d(i,j)=d(i,j)+(y0(i)-y0(j)+v*t(i,j)*(si(i)-si(j)))^2;endendendc=64-[d(2,1),d(3,1:2),d(4,1:3),d(5,1:4),d(6,1:5)];ceq=[];在MATLAB命令窗口计算如下：>>deitaini=zeros(1,6);>>vlb=-pi*ones(1,6)/6;vub=pi*ones(1,6)/6;>>OPT=optimset;OPT.LargeScale=’off’;>>[df,fval]=fmincon(@airfun,deltaini,[],[],[],[],vlb,vub,@airfunco,OPT);>>d1=dt*180/pi,fval=d1*d1结果为：d1=0.00000.00002.0683-0.4896-0.00551.5611fval=6.9548即各架飞机的调整角度分别为：0,0,2.0683，-0.4896，-0.0055,1.5611，调整角的平方和为6.9548。七、关于MATLAB的归纳与总结〔一〕关于MATLAB软件的总结在MATLAB下进行根本数学运算，只需将运算式直接打入提示号〔>>〕之后，并按入Enter键即可。MATLAB会将运算结果直接存入一变数ans，代表MATLAB运算后的答案〔Answer〕并显示其数值于荧幕上。">>"是MATLAB的提示符号，但在PC中文视窗系统下，由于编码方式不同，此提示符号常会消失不见，但这并不会影响到MATLAB的运算结果。MATLAB认识所有一般常用到的加〔+〕、减〔-〕、乘〔*〕、除〔/〕的数学运算符号，以及幂次运算〔^〕。另外，MATLAB将所有变数均存成double的形式。假设不想让MATLAB每次都显示运算结果，只需在运算式最后加上分号〔；〕即可，假设要显示变数y的值，直接键入y即可。另外，变数也可用来存放向量或矩阵，并进行各种运算。需要注意变数命名的规那么：第一个字母必须是英文字母；字母间不可留空格；最多只能有19个字母，MATLAB会忽略多余字母。MATLAB会忽略所有在百分比符号〔%〕之后的文字，因此百分比之后的文字均可视为程式的注解〔Comments〕。假设对MATLAB函数用法有疑问，可随时使用help来寻求线上支援helplinspace。对于循环来说，最简单的重复命令是for，其根本形式为：for变数=矩阵；运算式；end其中变数的值会被依次设定为矩阵的每一行，来执行介于for和end之间的运算式。因此,假设无意外情况，运算式执行的次数会等于矩阵的行数。最简单的逻辑命令是if,...,end，其根本形式为：if条件式；运算式；endifrand(1,1)>0.5,disp('Givenrandomnumberisgreaterthan0.5.');end假设要一次执行大量的MATLAB命令，可将这些命令存放于一个副档名为m的档案，并在MATLAB提示下键入此档案的主档名即可。此种包含MATLAB命令的档案都以m为副档名，因此通称M档案〔M-files〕。例如一个名为test.m的M档案〔所在的目录是d:\mlbook〕，包含一连串的MATLAB命令，那么只要直接键入test，即可执行其所包含的命令。严格来说，M档案可再细分为命令集〔Scripts〕及函数〔Functions〕。命令集的效用和将命令逐一输入完全一样，因此假设在命令集可以直接使用工作空间的变数，而且在命令集中设定的变数，也都在工作空间中看得到。函数那么需要用到输入引数〔Inputarguments〕和输出引数〔Outputarguments〕来传递资讯，这就像是C语言的函数,或是FORTRAN语言的副程序〔Subroutines〕。在前面提到的，test.m所在的目录是d:\mlbook。如果不先进入这个目录，MATLAB就找不到你要执行的M档案。如果希望MATLAB不管在何处都能执行test.m，那么就必须将d:\mlbook参加MATLAB的搜寻路径〔Searchpath〕上。要检视MATLAB的搜寻路径，键入path即可。此搜寻路径会依已安装的工具箱〔Toolboxes〕不同而有所不同。如果在每一次启动MATLAB后都要设定所需的搜寻路径，将是一件很麻烦的事。有两种方法，可以使MATLAB启动后，即可载入使用者定义的搜寻路径：1.MATLAB的预设搜寻路径是定义matlabrc.m〔在c:\matlab之下，或是其他安装MATLAB的主目录下〕，MATLAB每次启动后，即自动执行此档案。因此你可以直接修改matlabrc.m，以参加新的目录于搜寻路径之中。2.MATLAB在执行matlabrc.m时，同时也会在预设搜寻路径中寻找startup.m，假设此档案存在，那么执行其所含的命令。因此我们可将所有在MATLAB启动时必须执行的命令〔包含更改搜寻路径的命令〕，放在此档案中。每次MATLAB遇到一个命令〔例如test〕时，其处置程序为：将test视为使用者定义的变数。假设test不是使用者定义的变数，将其视为永久常数。假设test不是永久常数，检查其是否为目前工作目录下的M档案。假设不是，那么由搜寻路径寻找是否有test.m的档案。MATLAB储存变数的根本命令是save，在不加任何选项〔Options〕时，save会将变数以二进制〔Binary〕的方式储存至副档名为mat的档案，如下述：save：将工作空间的所有变数储存到名为matlab.mat的二进制档案。savefilename：将工作空间的所有变数储存到名为filename.mat的二进制档案。savefilenamexyz：将变数x、y、z储存到名为filename.mat的二进制档案。在数值微分方面，diff函数用以演算一函数的微分项，相关的函数语法有以下4个：diff(f)传回f对预设独立变数的一次微分值；diff(f,'t')传回f对独立变数t的一次微分值；diff(f,n)传回f对预设独立变数的n次微分值；diff(f,'t',n)传回f对独立变数t的n次微分值。数值微分函数也是用diff，因此这个函数是靠输入的引数决定是以数值或是符号微分，如果引数为向量那么执行数值微分，如果引数为符号表示式那么执行符号微分。在数值积分方面，int函数用以演算一函数的积分项，这个函数要找出一符号式F使得diff(F)=f。如果积分式的解析式(analyticalform,closedform)不存在的话或是MATLAB无法找到，那么int传回原输入的符号式。相关的函数语法有以下4个：int(f)传回f对预设独立变数的积分值；int(f,'t')传回f对独立变数t的积分值；int(f,a,b)传回f对预设独立变数的积分值，积分区间为[a,b]，a和b为数值式；int(f,'t',a,b)传回f对独立变数t的积分值，积分区间为[a,b]，a和b为数值式；int(f,'m','n')传回f对预设变数的积分值，积分区间为[m,n]，m和n为符号式。MATLAB解常微分方程式的语法是dsolve('equation','condition')，其中equation代表常微分方程式即y'=g(x,y)，且须以Dy代表一阶微分项y'D2y代表二阶微分项y''，condition那么为初始条件。对于求非线性方程式的实根，要求任一方程式的根有三步骤：先定义方程式。要注意必须将方程式安排成f(x)=0的形态，例如一方程式为sin(x)=3，那么该方程式应表示为f(x)=sin(x)-3。可以m-file定义方程式。代入适当范围的x,y(x)值，将该函数的分布图画出，借以了解该方程式的「长相」。由图中决定y(x)在何处附近(x0)与x轴相交，以fzero的语法fzero('function',x0)即可求出在x0附近的根，其中function是先前已定义的函数名称。如果从函数分布图看出根不只一个，那么须再代入另一个在根附近的x0，再求出下一个根。对于线性代数方程〔组〕求解，我们习惯将上组方程式以矩阵方式表示如下：AX=B，其中A为等式左边各方程式的系数项，X为欲求解的未知项，B代表等式右边之项。要解上述的联立方程式，我们可以利用矩阵左除\做运算，即是X=A\B。如果将原方程式改写成XA=B，其中A为等式左边各方程式的系数项，X为欲求解的未知项，B代表等式右边之项。注意上式的X,B已改写成列向量，A其实是前一个方