Matlab数据处理与分析1

/前言MATLAB一直是国际科学界应用和影响最广泛的软件工具，有着Mathematica和Maple无法比拟的优势和适用面。它不仅仅是一款数学软件，应用于微积分、概率统计、复变函数、线性变换、解方程、最优化、插值与数据显示等方面，也应用于模糊逻辑、小波分析、神经网络、图像处理、模式识别方面求解。另外数字信息处理、系统仿真、自动化、工程力学、信息与系统、模拟电路等方面都广泛的应用。使用好的数据处理方法和相应的软件工具对实验数据进行处理是大学理工数各专业学生应具备的基本技能，对于相关课程的学习也非常重要。MATLABT优点众多，本书只是重点讨论MATLAB在数据处理方面的应用。具体包括以下方面的内容：1．MATLAB安装与界面使用详解；2．数组与矩阵与其运算；3．MATLAB常用数值计算；4．代数方程与最优化；5．符号计算；6．插值与拟合；7．数据与函数的显示。在阅读过程中，要求上机执行书中的相关程序代码，熟练书的相关知识，要求勤查MATLAB自带的帮助系统，才能起到很好的学习效果。MATLAB的安装和界面使用本章详细讲述MATLABR2007a的安装和界面的使用。MATLAB安装首先双“setup.exe”之后，出现如图1-1的“WelcometotheMathWorksInstaller”窗口，然后点击“Next”按钮。1-1进入图1-2的“LicenseInformation”窗口，在Name框和Organization输入1-2相应名称（自行命名）。然后安装目录下的“serial.txt”文件，将序列号复制到“PleaseenteryourPersonalLicensePassword(PLP)”框中，然后点击“Next”按钮。待少许时间，防火墙会提示是否允许访问网络的提示框，如瑞星防火墙会出如下图1-3的对话框。1-3选择“总是允许”，并按“确定”。进入“MathWorksAccount”对话框如图1-4。图1-4选择“InstalltheproductversionsfrommyDVDorlocaldiskonly”选项按钮，进入“Licenseagreement”对话框，然后选择“Yes”按钮。图1-5进入“InstallationType”对话框，如图1-6所示，然后选择“Typical”按钮，点击“Next”。图1-6然后进入“”对话框，如图1-7所示。可以更改安装目录，也呆保持默认。如果系统盘的空间不空，则应更改安装目录，以保证机器的正确运行。图1-7在接下来的步骤中，皆选择“Next”按钮，即可。MATLAB的启动与界面详解在开始菜单中，选择：开始\程序\MATLAB\R2007a\matlabR2007a，如图1-8所示。图1-8然后进入MATLAB应用程序主界面，如图1-9所示。最上方是标题栏，显示版本信息：MATLAB7.4.0(R2007a)：再往下的一栏是菜单栏：图1-9接下来是工具栏：其中是指定当前工作路径，如果要向某文件夹中读取或保存相关内容，则通过点击，然后选择该文件夹作为当前工作路径。然后点击左侧的，则会出现CurrentDirectory对话框如图1-10所示，用来显示当前路径中的文件信息。它和Workspace处于同一具标签对话框中。图1-10右侧是CommandWindow窗口，如图1-11所示，是用来输入MATLAB指令的。一打开MATLAB则在CommandWindow中会出现“>>”，此符号表示MATLAB软件已准备好了，正在待命令的输入。如输入指令：A=[123]然后回车，则出现运算结果：A=123此指令表示输入一个一行三列的矩阵（即行向量），并将此矩阵（向量）保存在变量A中。关于矩阵与向量向量，在第二章中详细解释。图1-11指令执行完毕之后，会发现左下方的CommandHistory窗口中会多出现一行刚刚输入的指令。如图1-12所示，CommandHistory是用来保存输入过的命令，方便以后查找或再次使用。同时Workspace标签属性页也会多出一行。如图1-12所示，CurrentDirectory是用来显示可用的各种变量的。图中说是说明了刚才得到的变量A，它的值、最小、最大元素与类型等信息。图1-12图1-13在上任一处右击，然后在弹出菜单上选择需要查看的相关信息。如图1-14所示。如果需要用到以前的命令，可以在CommandHistory窗口中查到，可按日期查找，可以选择单条历史命令，也可以同时选择多条历史记录。然后右选择复制，将其粘贴到CommandWindow窗口中执行。也可将其作为文本复制到Word文档或txt文档中。如果双击历史记录中的命令，则系统会立即执行被双击的命令一次。如果CommandHistory窗口中一些命令不需要，可以一条或多条记录，然后按“Delete”键，将之删除。图1-14图1-15MATLAB的帮助详解点击主菜单中的Hellp\MATLABHelp或直接按F1，进入MATLAB自带的帮助界面，如图1-16所示。点击各个节点可以查看相关帮助信息。帮助文件里有详细的解释和丰富的实例。图1-16由帮助界面，可以看以MATLAB的强大功能，它可以与主流的软件开发工具进行混合编程开发。与Java和.Net无缝对接。只有勤查帮助或网上查找资料才能起到事半功倍的学习效果。数组与矩阵的概念与其运算在MATLAB中预定了一些符号，用于特定的含义，以下是一些预定义的符号：表2-1符号意义符号意义ans默认变量名pi圆周率eps机器可识别的最小的数flops浮点运算之数inf无穷大NaN非数i或j虚数单位realmax最大的实数realmin最小的实数intmax最大的整数其余的请参看联机帮助系统。数组与矩阵的概念MATLAB中，数组和矩阵本身是没有区别的，在内存中是一样的。只是针对不同的运算方式，将其为数组运算或矩阵运算。如果运算是按元素对应进行的，则称为数组运算。如果按线性代数学中的方式运算，则称为矩阵运算，如例2-1例2-1数组与矩阵的区别演示。设有矩阵A和B如下：以数组方式运算方式：，其中Inf表示无穷大（非数）MATLAB指令：A=[-11;02]%输入的矩阵保存在变量A中A=-1102B=[-11;02]%输入的矩阵保存在变量B中B=-1102A.*B%以数组方式相乘，ans=-1208A./B%以数组方式相除ans=-12Inf2几点说明：1、MATLAB中，矩阵（数组）按元素逐个输入的方法，就是将所有元素放在一对方括号内，行与行之间以分号“;”隔开，每一行中各元素之间以空格或逗号隔开。详细说明参看2.2节内容；2、MATLAB指令输入时，必需在英文状态下输入，否则会出错；3、如果一条指令以分号“;”作为结束符，则运算结果并不显示在CommandWindow窗口中，但保存在Workspace中。如果没有分号，则将结果显示在CommandWindow中；4、百分号“%”表示注释，从％开始到行未为方便人阅读所加的注释，不是指令中的一部分；5、注意数组方式运算有一个小黑点“.”；6、本书中凡MATLAB指令都用加粗字体；7、MATLAB中变量无需要声明就可直接使用，根据赋值符号“＝”右则表达式的类型来自动确定左侧的类型。如A=[-11;02]，A表示一个2阶矩阵。以矩阵方式运算（即按线性代数中的矩阵运算）：当然，矩阵的加法和减法两种方式都是一样的，MATLAB指令：A*B%矩阵方式相乘ans=-15-311A/B%矩阵方式相除ans=-1.00001.5000-3.00003.5000注意：如本书中后面内容只提与矩阵，只有当涉与到数组方式的运算时才将矩阵称为数组。矩阵的创建和操作在MATLAB中，有两种创建矩阵的方法。一是直接按元素逐个输入的方法，如例2-1所示；另一种就是使用MATLAB相关的指令来创建。一、直接输入元素创建直接按元素逐个输入来创建矩阵，就是将所有元素放在一对方括号内，行与行之间以分号“;”隔开，每一行中各元素之间以空格或逗号隔开。只有一行的矩阵称为行向量（也称为一维数组），只有一列的矩阵称为列向量。此方法可创建向量和矩阵。例2-2直接输入创建向量和矩阵演示。vr=[1234]%创建行向量，元素之间以空格隔开vr=1234vc=[1;2;3]%创建列向量，行之间以分号隔开vc=123m23=[123;456]%创建一个2×3行的矩阵m23=123456二、MATLAB指令创建在MATLAB中指令，更多时候也称为函数。可以使用MATLAB内置的函数来创建矩阵（数组）。以下以举例的方式说明。例2-3通过MATLAB指令创建向量和矩阵演示：1、指定起点：步长：终点。如果不指定步长，则将步长默认为1，最后一个元素不一定是终点，这取决于区间长度是否为步长的整数倍。该方法用于创建向量。v=0:0.2:1%以0为起点、1为终点、步长为0.2创建一个数组（行向量）v=00.20000.40000.60000.80001.0000v=0:pi%起点0、终点pi、默认步长1。最后一个元素不是终点。v=01232、linspace（起点，终点，元素个数），等分间隔。该方法用于创建向量。v=linspace(0,pi,3)v=01.57083.1416v=linspace(0,3,5)v=00.75001.50002.25003.00003．特殊矩阵的创建。创建特殊矩阵的常用函数：rand、magic、zeros、ones和eye等，需要深入研究请参看联机帮助。rand('state',0)%把均匀分布伪随机发生器置为0状态v=rand(2,3)%产生一个2×3的随机矩阵v=0.95010.60680.89130.23110.48600.7621m=magic(3)%产生一个3阶魔方矩阵m=816357492zeros(3)%产生一个3阶零矩阵ans=000000000zeros(2,3)%产生一个2×3的零矩阵ans=000000eye(2,3)%产生一个2×3的矩阵，左边2×2是一个单位矩阵ans=100010eye(3)%产生一个3阶单位矩阵ans=100010001ones(2,3)%产生一个元素全为1的2×3阶矩阵ans=111111另外还有其它特殊的矩阵创建函数，如有需要请参看帮助。矩阵的访问操作以下以举例的方式说明矩阵的访问操作如下：例2-4矩阵访问举例演示：v=[1234567];%生成一个行向量v(3)%查询第三个元素的值 ans=3v(3)=23%将第三个元素的值设为23v=12234567v([126])=[111216]%将下标为1、2、6的三元素的值设为11、12、16v=11122345167v(4:end)%查询第4至最后元素之间的所有元素ans=45167v(1:5)% 查询第1至5个元素ans=11122345m=[123;456]%产生一个新矩阵mm=123456m(2,3)%查询第2行第3列位置上的元素ans=6m(:,2)%查询第2列元素上所有行的元素ans=25m(2,:)%查询第2行上所有列的元素ans=456m(2,[12])%查询第2行上的第1、2列位置上的元素ans=45m(1,[23])=[8899]%将第1行上的第2、3列上的元素分别设为88和99m=18899456m(2,2)=518%将第2行2列位置上的元素设为518m=1889945186矩阵的基本运算操作矩阵的常用基本运算有加、减、乘、除、求逆等。数组方式和矩阵方式的运算符只差了一个小圆点，注意观察实例代码。1、加法和减法数组方式和矩阵方式都是一样的，就是直接将对应位置上的元素相加。如果是一个数和矩阵相加，则矩阵的每一个元素都加上这个数。例2-4矩阵加法操作演示：A=magic(3)%产生一个魔方矩阵A=816357492B=ones(3)%产生一个元素全为1的矩阵B=111111111A+B%矩阵的加法ans=92746851032+A%数与矩阵的加法ans=103857961142、乘法、除法与逆运算数组方式的乘法和除法是依元素对应相乘；矩阵方式的乘法则是按线性代数中的方法进行，矩阵的除是按线性代数中的取逆进行。左除：A/B，相当于A*B-1，右除：A\B，相当于A-1*B。逆运算按线性代数中的方法进行。例2-5矩阵乘、除操作演示A=magic(3);B=ones(3);A.*B%数组方式的乘法ans=816357492A*B%矩阵方式的乘法ans=151515151515151515M=[1,2;21]M=1221B=[1-1;10]B=1-110A./B%维数不匹配，将会出错???Errorusing==>rdivideMatrixdimensionsmustagree.M./B%数组方式的除法ans=1-22InfB./M%数据方式的除法ans=1.0000-0.50000.50000M/B%矩阵方式的除法ans=-23-13M*inv(B)%以矩阵方式运算，M乘以B的逆矩阵ans=-23-13inv(B)%计算B的逆ans=01-113、矩阵作为函数参数如果矩阵（数组）作为标准数学函数的参数，则对每一个元素都作同一函数计算。如V是一个行向量，R=sin(V)的运算结果R也是一个行向量，且R是的每一个元素都是由V中对应元素值求正弦值得到。例2-6矩阵作为函数参数演示v=[0pi/2pi3*pi/22*pi]%产生一个有4个元素的行向量v=01.57083.14164.71246.2832r=sin(v)%计算结果也是一个向量，和v的元素个相同r=01.00000.0000-1.0000-0.00004、常用的矩阵操作函数常用的矩阵操作函数如左右置换、上下置换、旋转，以下举例说明。例2-7矩阵常用操作演示A=magic(3)A=816357492flipud(A)%矩阵上下翻转ans=492357816fliplr(A)%矩阵左右翻转ans=618753294rot90(A)% 矩阵旋转90度ans=672159834非数、关系运算与逻辑操作1、非数NaN当表达式中如果分母出现零，或类似于的表达式运算时都会产生非数NaN，即Notanumber。NaN具有如下性质：NaN参与的运算结果也是NaN；非数没有大小的概念，因此不能将两非数去比较大小。非数真实的记录了运算结果，即数学中的无限变换趋，以下举例说明。例2-8非数使用演示log(0)%即相当于自变量趋向0时的极限为负无穷大ans=-Inft=1/0%相当于分母从右侧趋向趋向0时的极限为无穷大t=Infcos(t)%对非数计算余弦结果为非数ans=NaN非数要用于数据可视化中，如将图中某一指定部分镂空，将用到非数，详见第七章。2、关系运算与逻辑运算关系操作符有：==或eq(A,B),~=或ne,<或lt,>或gt,<=或le,>=或ge，以与&或and、|或or和~或nor等具体函数要求自行查阅帮助系统。另外有关函数如下：all,any,isqual,iempty,isfinite,isinf,isnan,isnumeric,isreal,isprime,isspace,isstr,ischar,isstudent,isunix,isvms,find.列表如下，详情请参看联机帮助。关系运算符如下表2-2关系运算符功能关系运算符功能<小于>=大于或等于<=小于或等于==等于>大于~=不等于例2-9关系运算与逻辑矩阵使用演示A=[159;347;268]A=159347268B=magic(3)B=816357492C=gt(A,B)％比较大小GreaterthanC=011000001whosC%查看C的详细信息NameSizeBytesClassAttributesC3x39logicalD=and(A,B)%求和运算D=111111111B>3%B中元素值大于1的位置对应1否则对应0，结果是一个逻辑矩阵ans=101011110B(find(B>3))%将B中元素值大于3的元素列出来ans=845967[r,c]=find(B>3)%元素值大于3的行号组成的数组r,列号组成数组cr=132312c=112233B.*(B>3)%B中不大于3位置上的元素设为零ans=806057490例2-10绘制[0，3pi]之间的曲线，并截去pi至2pi之间的曲线。x=linspace(0,3*pi);％自变量数组y=sin(x);％函数数组x1=(x<pi)|(x>2*pi);％逻辑数组y1=x1.*y;％截断数组plot(x,y1);％绘图小结、综合举例与练习例2-11综合举例，利用关系运算求近似极限，修补图形缺口。t=-2*pi:pi/10:2*pi;％自变量数组y=sin(t)./t;％函数值数组tt=t+(t==0)*eps;％修正后的自变量数组，元素值为零时，以最小机器数代替yy=sin(tt)./tt;％修正后的函数值数组subplot(1,2,1),plot(t,y),axis([-7,7,-0.5,1.2]),％绘制没有修正的图形subplot(1,2,2),plot(tt,yy),axis([-7,7,-0.5,1.2])％绘制修正后的图形图2-1矩阵的基本特征参数本章简要说明反映矩阵特征参数的一些量，如行列式、秩、条件数、范数、特征值与特征向量等问题。矩阵的基本参数以一下说明有矩阵信息的基本参数1、元素个数、行列数与其最大者、元素最大最小元素例3-1矩阵基本信息查询演示M=magic(3)M=816357492numel(A)%统计矩阵的元素个数ans=9size(M)%计算矩阵的行列数ans=33length(A)%计算行数与列数中的最大者ans=3max(M(:))%求出矩阵中所有元素中的最大者ans=9min(M(:))%求出矩阵中所有元素中的最小者ans=1矩阵的行列式、秩与范数计算行列式、秩与范数的指令分别是det、rank和norm例3-2矩阵行列式、秩与范数使用演示A=magic(3);det(A)%求A的行列式ans=-360rank(A)％计算矩阵的秩ans=3A=[569;351;861]A=569351861binf=norm(A,'inf')％计算无穷范数binf=20b2=norm(A,2)％计算2范数b2=15.4215条件数、矩阵的稳定性条件数是反映AX=b中，如果A或b发生细微变化，解变化的剧烈程度。如果条件数很大说明是病态方程方程，不稳定方程。例3-3矩阵条件数与稳定性演示A=[234;119;12-6]A=23411912-6con2=cond(A)％计算2-范式条件数con2=575.8240con1=condest(A)％计算1-范式条件数con1=817例3-4求解线性方程组：A=[234;119;12-6]％系数矩阵A=23411912-6b=[1;-7;9]％常数列b=1-79x=inv(A)*b％逆矩阵的方法求解x=11-1A\b％左除方法求解ans=11-1A=A+0.001％系数矩阵加上扰动A=2.00103.00104.00101.00101.00109.00101.00102.0010-5.9990b=b-0.001％常数列加上扰动b=0.9990-7.00108.9990x2=inv(A)*b％以逆矩阵的方法求解x2=0.95041.0297-0.9980特征值、特征向量与对角化与数学知识相关的概念请参考相关数学书籍例3-5特征值与特征向量演示：A=816357492E=eig(A)％计算特征值E=15.00004.8990-4.8990[V,D]=eig(A)％计算特征值组成的对角矩阵D和特征向量组成的矩阵VV=-0.5774-0.8131-0.3416-0.57740.4714-0.4714-0.57740.34160.8131D=15.00000004.8990000-4.8990注：有关正交化运算（orth函数）、三角分解（lu）、正交分解（qr）、特征值分解（eig）、奇异值分解（svd）的内容，请参看Matlab自还的帮助系统和相关数学书籍小结、综合举例与练习利用Matlab求解下列线性方程组：要求先分别列出系数矩阵和常数列向量，求出其行列式、秩、特征值与特征向量矩阵、2种条件数，系数矩阵的逆矩阵，然后求出其解，将系数矩阵再乘以解向量，看是否和常数向量相等。最后将系数矩阵加个0.0002的扰动，看其解是多少？微积分学中的基本求解问题符号运算简介符号运算以推理解析的方式进行，因此没有误差问题的困扰，能给出完全正确的解析解，如解析解不存在，则给按任意指定的精度给出数值解。符号运算使用起来非常简单，与大学数学中的书写类似。符号函数的定义使用指令sym与syms，以下代码段举例说明其应用。例4-1sym与syms的使用y=sym('2*sin(x)*cos(x)')％通过字符串来创建符号对象，自变量为xy=2*sin(x)*cos(x)y=simple(y)％将符号函数对像化简y=sin(2*x)symsab;％声明符号变量a,by=sin(a)*cos(b)-cos(a)*sin(b);％符号对象yy=simple(y)％将符号对象y化简y=sin(a-b)A=sym('[12;32]')％通过字符串创建矩阵符号对象AA=[1,2][3,2]da=det(A)％计算A的行列式da=-4ia=inv(A)％计算A的逆矩阵，注意结果含有分数ia=[-1/2,1/2][3/4,-1/4]ea=eig(A)％计算A的特征值ea=4-1[ev,ea]=eig(A)％计算特征值和特征向量。ev=[-1,1][1,3/2]ea=[-1,0][0,4]注意：使用符号运算就像平时数学运算一样方便。另外MATLAB中还有针对多项式操作的一系列函数：colllect：同幂次的项系数进行合并；expand：按多项式、三角函数或指数对数函数等展开；factor：因式分解；horner：将多项式分解成嵌套形式；simple：将表达式化简。如需要时可查看帮助，此不再详细展开。符号对象的精度控制在符号计算中，当符号常数或符号结果需要以数值形式给出时，可以灵活地按指定精度输出数值。与精度有关的三个函数指令：double，digits，vpadigits：显示当前采用的数值计算精度；digits(n)：设置今后数值计算以n位相对精度进行；xs=vpa(x)：在digits指定精度下，给出x的数值型结果；xs=vpa(x,n)：在n位相对精度下，给出x的数值型结果。以下举例说明。例4-2精度控制演示。vpa(pi)％依当前精度输出ans=3.1415927eval('pi')％依默认精度输出ans=3.1416vpa(pi,'100')％依指定精度输出ans=3.149323846264338327953993752117068a=sym('1/3')％创建符号对象a=1/3digits(2);％设置当前精度为2vpa(a)％以当前精度输出ans=.33digits(19);％设置当前精度为19vpa(a)％以当前精度输出ans=.3333333333333333333极限求解MATLAB中计算极限的函数是limit，其语法如下：limit(F,x,a)：计算limit(F,x,a,'right')：limit(F,x,a,'left')：例4-3求符号极限演示：symsxath％说明字符变量 limit(sin(x)/x)％计算极限，默认求右趋于零的极限ans=1limit(1/x,x,0,'right')％对x求右趋于零的极限ans=Inflimit(1/x,x,0,'left')％对x求左趋于零的极限ans=-Inflimit((sin(x+h)-sin(x))/h,h,0)％对h求趋于零的极限ans=cos(x)v=[(1+a/x)^x,exp(-x)]％v=[(1+a/x)^x,exp(-x)]limit(v,x,inf,'left')％对x求左趋于正无穷大的极限ans=[exp(a),0]导数与其几何应用导数问题分两类，一类是已知函数解析式，然后计算其导数解析式或计算某点的导数值，此类问题很好处理。另一类问题是只有实验得到的离散数据，要得到各阶导数值，该类问题比较麻烦，常用的方法是中心差分法。几何上介绍点积、叉积、混合积，切平面和法线、梯度场，流线场。一、解析法求导数问题解析法求导使用符号运算，求导函数为diff，语法如下：diff(S);％对S求导，依26字母中最接近x的字母为自变量diff(S,'v');％以v为自变量，对S求导diff(S,'v',n);％以v为自变量，对S求n阶导数例4-3求符号极限演示：da=diff('x^5+3*x+5')％求导数da=5*x^4+3f=sym('log(x)/exp(x^2)')％定义符号函数 f=log(x)/exp(x^2)simplify(diff(f))％对函数求导数后再化简ans=-exp(-x^2)*(-1+2*log(x)*x^2)/xsimplify(diff(f,2))％对函数求二阶导数后化简ans=exp(-x^2)*(-1-4*x^2+4*log(x)*x^4-2*log(x)*x^2)/x^2二、向量的点积、叉积和混合积点积、叉积与混合积的相关概念，请参见数据书籍。例4-4点积、叉积与混合演示：x1=[1980]x1=1980x2=[2697]x2=2697y=dot(x1,x2)％点积运算y=128y1=cross(x1,x2)％维数是3才能进行叉积运算???Errorusing==>crossat37AandBmusthaveatleastonedimensionoflength3.x1=[952]x1=952x2=[327]x2=327xd=dot(x1,x2)％计算点积xd=51xcr=cross(x1,x2)％计算叉积xcr=31-573xcr2=cross(x2,x1)％计算叉积xcr2=-3157-3x3=[247]x3=247ydc=dot(x3,cross(x1,x2))％混合积运算，用于计算平行六面体的体积ydc=-145各种积分问题积分有符号方法和数值方法，以下分别说明。符号积分的函数是int，它的语法是：int(S);％计算S的一个原函数。int(S,v);％以v为积分变量计算S的一个原函数int(S,a,b);％计算符号积分，a与b是下限和上限，S是被积函数int(S,v,a,b);％计算符号积分，a与b是下、上限，S是被积函数，v是积分变量例4-5符号积分演示：f=sym('sin(s+2*x)')％定义符号函数f=sin(s+2*x)int(f)％求原函数ans=-1/2*cos(s+2*x)int(f,'s')％以s为积分变量求原函数ans=-cos(s+2*x)int(f,pi/2,pi)％计算指定区间的定积分ans=-cos(s)int(f,'s',pi/2,pi)％指定积分变量，指定区间求定积分ans=-sin(2*x)+cos(2*x)int(f,'a','b')％指定区间上求定积分ans=1/2*cos(s+2*a)-1/2*cos(s+2*b)例4-6求与。F1=int('1/log(t)','t',0,'x')％计算符号积分，无初等解析式。Warning:Explicitintegralcouldnotbefound.>Inat58Inat9F1=int(1/log(t),t=0..x)F1=int('1/log(t)','t',0,'1/2')％计算积分，无初等解析式。F1=-Ei(1,log(2))vpa(F1)％以指定精度得到数值解。ans=-.378671672720718463656例4-7求积分。注意：内积分上下限都是函数。symsxyzF2=int(int(int(x^2+y^2+z^2,z,sqrt(x*y),x^2*y),y,sqrt(x),x^2),x,1,2)VF2=vpa(F2)F2=1610027357/656/348075*2^(1/2)+14912/4641*2^(1/4)+64/225*2^(3/4)VF2=224.9232805以下实例是数值计算方法计算积分。例2-8绘制下列空间曲线的图形并计算长度，参数方程：，长度：t=0:0.05:3*pi％生成数组plot3(sin(2*t),cos(t),t)%绘制三维曲线图4-1f=inline('sqrt(4*cos(2*t).^2+sin(t).^2+1)')%d定义一个内联函数f=Inlinefunction:f(t)=sqrt(4*cos(2*t).^2+sin(t).^2+1)%functionlen=quad(f,0,3*pi)％使用数值方法计算定积分的值len=17.2220级数问题求解级数求各的函数是symsum，它的语法形式是：r=symsum(s,a,b);r=symsum(s,v,a,b);以下举例说明它们的使用方法，详细请查阅MATLAB自带的帮助文件。例4-7求和与，数值转换与精度控制演示：A=symsum(sym('(2*n-1)^2'),1,'N')％求和A=11/3*N+8/3-4*(N+1)^2+4/3*(N+1)^3simplify(A)％化简ans=-1/3*N+4/3*N^3B=symsum(sym('1/((2*n-1)^2)'),1,inf)％求无穷项和B=1/8*pi^2％digits(15)％设置当前精度vpa(B)％必须使用vpa函数才能用相应精度输出ans=1.23370055013617例4-10求和与和，数值转换与精度控制演示：symsxkr=symsum(1/k^2,1,20)％1至20求和r=17299975731542641/170720r=symsum(1/k^2,10,20)％10至20求和r=32561/54192375991353600代数方程与最优化代数方程求解代数方程的函数是solve，它的语法是：solve(eq)solve(eq,var)solve(eq1,eq2,...,eqn)g=solve(eq1,eq2,...,eqn,var1,var2,...,varn)说明：输入参数符号表达式或字符串。如果输入参数是形如的字符串，即：solve()，则执行指令求解方程。可以求解单个方程，也可求解方程组。参数var，var1，var2，...，varn用于指定未知变量，eq1，eq2，...，eqn是用符号表达式或字符串表示的方程。另外MATLAB还提供了函数fsolve，可用于求解多元方程的一个实根，详细使用方法请参看帮助文件，此处主要介绍solve.例5-1求解方程，其中分别将x和b视为未知数。A=solve('a*x^2+b*x+c')％求解方程，A表示解组成的数组，x为未知数。A=-1/2*(b-(b^2-4*a*c)^(1/2))/a-1/2*(b+(b^2-4*a*c)^(1/2))/aA=solve('a*x^2+b*x+c','b')％求解方程，b为未知数。A=-(a*x^2+c)/x例5-2求解方程组，和S=solve('x+y=1','x-11*y=5')％求解方程组S=x:[1x1sym]y:[1x1sym]S.x,S.y％显示结果ans=4/3ans=-1/3A=solve('a*u^2+v^2','u-v=1','a^2-5*a+6')％解方程组并显示结果的结构A=a:[4x1sym]u:[4x1sym]v:[4x1sym]A.a,A.u,A.v％显示结果的A的各个部分ans=2233ans=1/3+1/3*i*2^(1/2)1/3-1/3*i*2^(1/2)1/4+1/4*i*3^(1/2)1/4-1/4*i*3^(1/2)ans=-2/3+1/3*i*2^(1/2)-2/3-1/3*i*2^(1/2)-3/4+1/4*i*3^(1/2)-3/4-1/4*i*3^(1/2)a=A.a;u=A.u;v=A.v;％显示结果的A的各个部分分别保存在a,u,v中[eval('a.*u.^2+v.^2')eval('u-v-1')eval('a.^2-5.*a+6')]％检验结果是否正确ans=[2*(1/3+1/3*i*2^(1/2))^2+(-2/3+1/3*i*2^(1/2))^2,0,0][2*(1/3-1/3*i*2^(1/2))^2+(-2/3-1/3*i*2^(1/2))^2,0,0][3*(1/4+1/4*i*3^(1/2))^2+(-3/4+1/4*i*3^(1/2))^2,0,0][3*(1/4-1/4*i*3^(1/2))^2+(-3/4-1/4*i*3^(1/2))^2,0,0]a.*u.^2+v.^2％检验结果是否为零，下面的结果似乎表示方程左边不为零ans=2*(1/3+1/3*i*2^(1/2))^2+(-2/3+1/3*i*2^(1/2))^22*(1/3-1/3*i*2^(1/2))^2+(-2/3-1/3*i*2^(1/2))^23*(1/4+1/4*i*3^(1/2))^2+(-3/4+1/4*i*3^(1/2))^23*(1/4-1/4*i*3^(1/2))^2+(-3/4-1/4*i*3^(1/2))^2u-v-1％检验结果是否为零，结果显示方程左边为零ans=0000a.^2-5.*a+6％检验结果是否为零，结果显示方程左边为零ans=0000eval(a.*u.^2+v.^2)％将符号变量转成数值，结果显示方程左边几乎可看成零ans=1.0e-015*000.11100.1110无条件最优化问题求解MATLAB中解决无条件极值的函数有两个：fminunc和fminsearch，其中fminunc算法效率要优于fminsearch,此处主要讲解fminunc。它主要是解决最小值问题（最大值问题也可化为最小值问题）：，其中x是向量，f(x)是目标函数。它的语法形式如下：x=fminunc(fun,x0)x=fminunc(fun,x0,options)[x,fval]=fminunc(...)函数使用说明：1、fminunc尝试找到标量函数fun的最小值，从x0开始始搜索。x0是预先指定的初值，fun是目标函数，返回值x即为所搜索到的局部最小值点，fval是函数fun在局部最小值点x处的函数值。还有其它返回值，具体请参看帮助文件。2、该函数搜索到的可能是局部极小值点，而不是全局极小值。很难找到万能的搜索全局最小值的算法函数，因此在使用该函数可考虑由函数图像中选择一个极小值点附近处的某值x0作为初值，才能确定全局最小值。3、函数fun可以是内联形式的函数，也可以是以M文件形式存在的函数，还可以是无名函数，因内联形式和无名函数形式相对简单，本节只举例讲述此两者，至于M文件形式的函数，请参看帮助文件。4、有时搜索的速度较慢，则需要使用函数梯度，可加速搜索。但略复杂，如实际需要可参看帮助文件。5、options的关键参数有5个分别是：options.Display：用于显示设置，有四个选取值off、iter、final和notify；options.MaxFunEvals：允许函数计算的最多次数；options.MaxIter:允许的最多迭代次数；options.TolFun:函数值计算的终止容差，缺少值1.0000e-004；options.TolX：自变量计算的终止容差，缺少值同上。以后涉与options参数的同上。例5-3求函数的最小值以下指令是直接求解，不显示计算的迭代过程。f=inline('3*x(1)^2+2*x(1)*x(2)+x(2)^2');％内联形式的函数，注意形式x0=[1,1];％初值[x,fval]=fminunc(f,x0)x=1.0e-006*0.2541-0.2029fval=1.3173e-013以下代码是显示中间的迭代过程f=inline('3*x(1)^2+2*x(1)*x(2)+x(2)^2');x0=[1,1];fopt=optimset;％设置选项，用于显示执行过程fopt.Display='iter';％显示迭代过程[x,fval]=fminunc(f,x0,fopt)IterationFunc-countf(x)Step-sizeoptimality0368160.250.1251290.080576310.3893120.051033410.2764150.0023402910.1585180.00011471610.03666213.4414e-00710.001167249.12732e-01014.45e-0058271.31734e-01311.16e-006x=1.0e-006*0.2541-0.2029fval=1.3173e-013例5-4求函数的最小值，与f=@(x)sin(x)+3;％无名函数形式[x,f]=fminunc(f,4)x=4.7124f=2.0000f=inline('sin(x(1))+cos(x(2))+3');％内联形式，x0=[45];[x,fv]=fminunc(f,x0)x=4.71243.1416fv=1.0000线性规划求解MATLAB求解线性规划的函数是linprog，它是计算如下线性规划问题的最小值：，其中注意，等式条件(2)中，是矩阵，是向量，不等式条件(1)中，A是矩阵，b是向量，不等式条件(3)中的lb和ub也是向量。linprog的语法是：x=linprog(f,A,b)x=linprog(f,A,b,Aeq,beq)x=linprog(f,A,b,Aeq,beq,lb,ub)x=linprog(f,A,b,Aeq,beq,lb,ub,x0)x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options)说明：f即是目标函数，A,b,Aeq,beq,lb,ub,x0同前所述，x是最小值点，fval是最小值点处的函数值（即搜索到的最小值）。如果没有相应的条件，则对应的矩阵或向量应设为空矩阵（空数组），以下举例说明。例5-5求函数的最小值，其中x满足条件：(1)先按变量的顺序排好，然后用系数表示目标函数，如此处：f=[-5;-4;-6];(2)由于没有等式条件，所以Aeq、beq都是空矩阵即[]；(3)不等式条件中的第三个式子中，没有x3，所以系数为零，故有：(4),,由于没有上限要求，故设为无穷大。解：f=[-5;-4;-6];％目标函数的系数A=[1-11324320];％不等式各变量系数b=[20;42;30];％不等式条件中的常列lb=[0;0;0];％各变量的下限ub=[inf;inf;inf];％各变量的上限[x,fval]=linprog(f,A,b,[],[],lb,[]);％求解运算Optimizationterminated.x,fval％查看结果x=0.000015.00003.0000fval=-78.0000注意：凡一如果是大于的不等式，则将其转化成小于不等式，求最大值的，转化成求最小值。练习1：求解函数的最大值，其中x满足条件：f=[-3/4,150,-1/50,6];％目录函数的系数，取反是为了得到最小值，即原函数的最大值参考代码：Aeq=[];beq=[];％没有等式条件，所以相应的矩阵为空A=[1/4,-60,-1/50,9;1/2,-90,-1/50,3];％不等式条件中的各变量系数矩阵b=[0,0];％不等式条件中的常数列lb=[-5;-5;-5;-5];ub=[inf;inf;1;inf];％各自变量的下限和上限x0=[1;1;1;1];％初始值[x,fval]=linprog(f,A,b,Aeq,beq,lb,ub,x0);％求解运算Optimizationterminated.x,fvalx=-5.0000-0.19471.0000-5.0000fval=-55.4700说明：1、对于目标函数是二次型的形式且约束条件仍然为线性的，则应使用专门的相应二次型规划函数：quadprog；2、对于有约束的非线性最优化的问题，应使用非线性规划问题求解函数：fmincon.3、对于0-1规划，即自变量只能取0或1两个值的规划问题。MATLAB中相应的规划函数是bintprog。由于篇幅原因，此两类问题不再详述，如有需要请参看帮助文件。数据与函数的可视化显示视觉是人们感受世界、认识自然的最重要依靠。数据可视化的目的在于：通过图形，从一堆杂乱的离散数据中观察数据间的内在关系，感受由图形所传递的内在本质。MATLAB一向注重数据的图形表示，并不断地采用新技术改进和完备其可视化功能。本章将系统地阐述：离散数据表示成图形的基本机理；曲线、曲面绘制的基本技法和指令；重点与关键步骤：一、准备数据二、显示数据二维显示常用的二维绘图函数有plot、quiver、contour。除此之外还有不太常用的绘图函数，如polar、loglog、semilogx、semilogy、bar、barh、rose、stairs、area、pie、scatter、hist、erroerbar、stem、feather、commet、compass等，详情请看参看帮助文件这里只介绍plot、quiver和contour。一、plot绘图函数简介plot的语法如下plot(Y)plot(X1,Y1,...)plot(X1,Y1,LineSpec,...)plot(...,'PropertyName',PropertyValue,...)plot(axes_handle,...)h=plot(...)函数说明：1、当输入参数只有一个时，如plot(Y)，Y是输入向量，以Y的下标作为横坐标，以Y的值作为纵坐标绘图，默认的方式是各离散点之间使用线段连接起来，当点取得足够密时，看起来就是一条光滑的曲线。2、plot(X1,Y1,...,Xn,Yn)针对每一对坐标数据参数(Xi,Yi)绘制一条曲线，Xi,Yi都一向量。每一对坐标数据参数后可紧跟一个控制线型或颜色的的参数，也可没有。如果没有线型或颜色的参数，则采用默认的线型和颜色绘制曲线。注意，Xi,Yi也可以是矩阵，情况略复杂，成对的矩阵参数必需是同型矩阵（行列数都相同），以Xi的某列为横坐标向量，同Yi中相同列号的列为纵坐标向量绘制曲线，具体请看实例演示。例6-1绘制函数在上的不同参数下的图像代码段1：x=-pi:pi/10:pi;％第一步准备数据：自变量数组y=tan(sin(x))-sin(tan(x));％第一步准备数据（函数数组）plot(x,y,'--rs','LineWidth',2,...％一条指令要多行显示，则在行后加三个圆点“…”'MarkerEdgeColor','k',...％各点连线的颜色，'MarkerFaceColor','g',...％坐标点的颜色'MarkerSize',10)％'--rs'中的--表示连线线型，rs表示点形状得到如图6-1所示的图像图6-1图6-2代码段2x=-pi:pi/10:pi;y=tan(sin(x))-sin(tan(x));plot(x,y,'--','LineWidth',2,...％和代码段相比，只是坐标点的形状变了。'MarkerEdgeColor','k',...'MarkerFaceColor','g',...'MarkerSize',10)此代码段产生图6-2的图像。注：b蓝,g绿,r红,c青,m:品红,y:黄,k:黑,w:白例6-2绘制函数在上的图像图6-3图6-4x=-pi:.1:pi;y=sin(x);plot(x,y);％第二步绘图，参数分别对应自变量和函数set(gca,'XTick',-pi:pi/2:pi)％设置x轴的刻度区间长度set(gca,'XTickLabel',{'-pi','-pi/2','0','pi/2','pi'})％标示刻度例6-3二维曲线绘制，并显示网格：t=0:pi/100:2*pi;y=sin(t);％第一步准备数据plot(t,y)％第二步绘图gridon%Turnongridlinesforthisplot得到图6-4的图像。例6-3二维曲线平移正弦曲线组绘制：t=0:pi/100:2*pi;y=sin(t);y2=sin(t-0.25);％向右平衡0.25单位，绘图区间不变y3=sin(t-0.5);％向右平衡0.5单位，绘图区间不变plot(t,y,t,y2,t,y3);holdon;％保留之前所画图，之后的画图时不清除以前画的图plot(t,y-1,t,y2+1,t,y3);％上、下各平衡1个单位得到图6-5，关于holdon的详情，请参看帮助文件图6-5图6-6例6-4一组椭圆的绘制演示：。th=[0:pi/50:2*pi]';％参数值数组101×1a=[0.5:.5:4.5];％轴长数组1×9X=cos(th)*a;％水平坐标值数组101×9Y=sin(th)*sqrt(25-a.^2);％纵坐标值数组101×9plot(X,Y)％绘图，以X,Y中对应的列分别的曲线的横、纵坐标绘曲线。因X,Y共9列，所以共绘制9条曲线。得到图6-6例6-5一组正弦曲线的绘制演示：t=(0:pi/50:2*pi)';％自变量数据准备101×1，注意“'”表示取转置k=0.4:0.1:1;％缩放系数数组1×7Y=cos(t)*k;％函数的相关数据101×7plot(t,Y)％绘图注：相当于：plot(t,0.5*cos(t),t,0.6*cos(t),t,0.7*cos(t),t,0.8*cos(t),t,0.9*cos(t),t,cos(t))图6-7图6-8例6-6用图形表示连续调制波形与其包络线演示：t=(0:pi/100:pi)'; %101×1y1=sin(t)*[1,-1]; %101×2，包络y2=sin(t).*sin(9*t); %101×1，曲线t3=pi*(0:9)/9; %1×10，零点坐标y3=sin(t3).*sin(9*t3);％零点函数值，1×10plot(t,y1,'r:',t,y2,'b',t3,y3,'bo') axis([0,pi,-1,1]) %坐标轴比例大小注：axis([xminxmaxyminymax])：setsthelimitsforthex-andy-axisofthecurrentaxes.练习：绘制有的平移曲线组，要求有五条曲线练习二：使用Matlab在一个plot指令中绘制y=sinx,y=0.2sinx,y=0.4sinx,y=0.6sinx,y=0.8sinx的组图，大致如下图6-9所示：图6-9图6-10二、contour绘图函数简介contour用于绘制三维曲面的二维等高线，并带有标注。它的常用的基本语法如下：contour(X,Y,Z)contour(X,Y,Z,n)contour(X,Y,Z,v)函数说明：X,Y,Z是同型矩阵（行列数相同），分别表示记录空间曲面网格点的,x,y,z坐标矩阵，n表示等高线的条数，是一个标量。v是一个向量，其元素值表示高度值（即z坐标值），在v指定的z坐标高度画等高线。当然还有其它语法形式，具体参看帮助文件。由于空间曲面的描述需要网格点，而产生网格数据的函数是meshgrid，具体看帮助系统，此处仅举例说明。例6-7绘制曲面的等高线，其中-2<x<2,-2<y<2。[X,Y]=meshgrid(-2:.2:2,-2:.2:3);%X、Y分别是记录x和y坐标的矩阵Z=X.*exp(-X.^2-Y.^2);%Z分别是记录z坐标的矩阵[C,h]=contour(X,Y,Z);%绘制等高线，C记录各等高线数据的矩阵，h是等高线图形对象colormapcool%再改颜色风格set(h,'ShowText','on');%显示条等高线的标注注意：可以在Workspace中查看C,h，上述代码得到图6-10。三、gradient与quiver绘制二维梯度场1、gradient用于绘制二元函数的梯度场，其运算是基于离散网格数据点的，梯度使用差分来近似表示梯度的理论值。其常用的语法形式如下：FX=gradient(F)[FX,FY]=gradient(F)[FX,FY,FZ,...]=gradient(F)[...]=gradient(F,a,b)函数说明：F表示二元函数的函数值矩阵，即Z=f(x,y)中的z坐标值的矩阵。FX是水平方向的偏导数值矩阵，即水平方向的差分，FY表示垂直方向的差分值矩阵。如是三元函数，则FZ表示z轴方向的差分值矩阵。a,b用于指定计算差分时两个方向的步长。如果参数只两个，一个是F,另一个是标量参数，则两个方向步长相同。2、quiver函数是根据gradient函数得到的水平和垂直梯度值来绘制梯度场，有时为了方便，将梯度场和等高线绘制在同一幅图中。quiver的常用语法如下：quiver(x,y,u,v)h=quiver(...)x,y,u与v是同型矩阵，x与y记录定义域中网格点的横坐标值与纵坐标值，u与v记录网格点的水平与纵向梯度值，h表示所绘梯度场对象。例6-8绘制曲面的梯度场与带有等高线的梯度场，其中-2<x<2,-2<y<2。x=-2:.2:2;%定义域中x方向的定义域y=-1:.2:1;%定义域中y方向的定义域[xx,yy]=meshgrid(x,y);%产生网格矩阵以记录网格点的x坐标和y坐标的值zz=xx.*exp(-xx.^2-yy.^2);%记录网格点的z坐标[px,py]=gradient(zz,.2,.2);%得到x,y方向梯度quiver(x,y,px,py,2);%绘制梯度场，第5个参数表示箭头的长度。此段代码得到图6-11图6-11图6-12以下是代码段产生一个带有等高线的梯度场，如图6-12所示。v=-2:0.2:2;[x,y]=meshgri