[工学]matlab教案
54目录第1章绪论21.1科学计算简介21.2Matlab概述2第2章Matlab运算基础52.1变量与赋值52.2矩阵62.3表达式92.4数学函数9第3章Matlab程序设计103.1m文件103.2数据的输入输出103.3程序结构113.4函数文件12第4章图形与声音144.1二维图形144.2三维图形144.3图形窗口控制154.4图形控制164.5动画164.6声音16第5章线性代数185.1矩阵185.2向量空间195.3线性方程组215.4特征值与特征向量255.5二次型及其标准型*27第6章数据处理与多项式306.1基本统计处理306.2多项式316.3数据插值326.4曲线拟合356.5离散傅立叶变换36第7章数值积分与微分方程387.1数值积分387.2数值微分387.3常微分方程的数值解417.4非线性方程(组)求解457.5函数优化47第8章符号计算488.1符号计算基础488.2微积分498.3线性代数518.4方程求解53\n54第1章绪论1.1科学计算简介科学计算，即对科学和工程中的数学问题进行数值计算。数值计算的过程主要包括建立数学模型、建立求解的计算方法、计算机实现三个阶段。数值计算的特点是计算方法比较复杂，方法种类多种多样，如数值微分、数值积分、常/偏微分方程、线性代数方程、有限元等。数值计算所关心的焦点是计算精度（误差影响）。科学计算可分为两类：一类是纯数值的计算，例如求函数的值，方程的数值解;另一类计算是符号计算，又称代数运算，这是一种智能化的计算，处理的是符号。符号可以代表整数，有理数，实数和复数，也可以代表多项式，函数，还可以代表数学结构如集合，群的表示等等。我们在数学的教学和研究中用笔和纸进行的数学运算多为符号运算。主要的数学软件有：Mathematica、MATLAB、Maple和MathCAD。尽管计算机代数系统在代替人进行繁琐的符号运算上有着无比的优越性，但计算机毕竟是机器，它只能执行人们给它的指令。此外，虽然计算机代数系统包含大量的数学知识，但这仅仅是数学的一小部分，目前有许多数学领域计算机代数系统还未能涉及。1.2Matlab概述1980年前后，当时的新墨西哥大学计算机系主任CleveMoler教授在讲授线性代数课程时，发现了用其他高级语言编程极为不便，便构思并开发了MATLAB(MATrixLABoratory，即矩阵实验室),这一软件利用了当时数值线性代数领域最高水平的EISPACK和LINPACK两大软件包中可靠的子程序，用Fortran语言编写了集命令翻译、科学计算于一身的一套交互式软件系统。只要给出一条命令，立即就可以得出该命令的结果，而无需像C和Fortran语言那样，首先编写源程序，然后对之进行编译、连接，最终形成可执行文件。这无疑会给使用者带来了极大的方便。此后，CleveMoler等人成立了MathWorks的公司，用C语言作了完全的改写，于1984年推出了第一个MATLAB的商业版本。其后又增添了丰富多彩的图形图像处理、多媒体功能、符号运算和它与其他流行软件的接口功能，使得MATLAB的功能越来越强大，应用范围越来越广。到目前为止其最高版本7.0版已经推出。MATLAB具有强大的数学运算能力、方便实用的绘图功能及语言的高度集成性，它在其他科学与工程领域的应用也是越来越广，并且有着更广阔的应用前景和无穷无尽的潜能。MATLAB可以将使用者从繁琐、无谓的底层编程中解放出来，把有限的宝贵时间更多地花在解决问题中，这样无疑会提高工作效率。目前，MATLAB已经成为国际上最流行的科学与工程计算的软件工具，现在的MATLAB已经不仅仅是一个“矩阵实验室”了，它已经成为了一种具有广泛应用前景的全新的计算机高级编程语言了，在国内外高校和研究部门正扮演着重要的角色。可以预见，在科学运算、自动控制与科学绘图领域MATLAB语言将长期保持其独一无二的地位。一、MATLAB的主要特点l有高性能数值计算的高级算法特别适合矩阵代数领域；l有大量事先定义的数学函数并且有很强的用户自定义函数的能力；l有强大的绘图功能以及具有教育、科学和艺术学的图解和可视化的二维、三维图；l适合个人应用的强有力的面向矩阵(向量)的高级程序设计语言；\n54l与其它语言编写的程序结合和输入输出格式化数据的能力；l有在多个应用领域解决难题的工具箱。本教程基于MATLAB6.5版。一、操作环境matlab主界面：菜单、工具栏、工作区（命令窗口）、workplacebrower、commmandhistory等。Figtrue窗口：m文件编辑窗口：有些工具箱也提供了操作界面，可以通过命令或选择start的toolboxes菜单进入。二、引例>>2+6–4ans=4>>ans/2ans=2也可以定义变量>>a=5a=5>>b=6b=6\n54>>c=b/ac=1.2000可以使用内建的变量和函数>>pians=3.1416>>sin(ans/4)ans=0.7071一、帮助在MATLAB系统中相关的线上(on-line)求助方式有：1.利用help或helpwin指令：在命令窗口键入help
，例如helpsqrt,help('sqrt'),helpwinsqrt。2.利用lookfor指令：键入的lookfor(key-word)，列出所有相关的题材，例如lookforcosine,lookforsine。3.利用指令视窗的功能选单中的Help，从中选取TableofContents（目录）或是Index（索引）。4.利用doc指令：doc；或docfunction\n54Matlab运算基础1.1变量与赋值一、变量matlab变量主要是矩阵，矩阵的元素可以是数值或字符串。其中数值可以是实数或复数，虚数单位用i或j表示。如a=3+9i（最好不要写成3+9*i，因为i、j是可以作为普通变量被赋值的）特别地，如果矩阵只有一行或一列，就称为向量(vector),，如果只有一行一列，就称为标量(scalar)。其他的变量类型有：多维数组、结构、对象等。变量的命名：matlab变量无需声明即可直接使用，变量名以字母开头，接字母、数字或下划线，最多63个字符（多余的被忽略），区分大小写。系统变量：ans：valueofanexpressionwhenthatexpressionisnotassignedtoavariableeps：floatingpointprecision2.2204e-016(2-52)pi：π,(3.141492...)realmax：最大正浮点数1.7977e+308(21024)realmin：最小正浮点数2.2251e-308(2-1022))Inf：∞，正无穷大,anumberlargerthanrealmax,theresultofevaluating1/0.NaN：notanumber,theresultofevaluating0/0i,j：虚数单位不要对系统变量赋值，否则使用内建函数时会出错。例外：i,j常作为循环变量使用。二、赋值语句1）变量=表达式2）直接输入表达式：结果被赋给系统变量ans注意：l语句后加分号,表示仅执行，不显示结果l续行符为…（三个句点）l注释：以%开头例子：p14，计算表达式的值，并将结果赋给变量a。>>a=2*sin(85*pi/180)/(1+sqrt(5)+3i)a=0.3311-0.3070i变量的查看：who,whos。who列出工作空间中的所有变量名，whos同时给出变量的维数和性质。变量也列在workplacebrowser中。数据的格式：format格式符，格式符参见p18。如formatshort。也可以通过file->prefrences菜单，在prefrences对话框的commandwindow项中设置。\n54变量的清除：clear，清除工作空间中的所有变量，等同edit->clearworkplace菜单。清除命令窗口：clc，等同edit->clearconmmandwindow菜单。1.1矩阵一、矩阵的建立1）直接输入将矩阵元素用方括号括起来，同一行各元素间用空格或逗号隔开；各行之间用分号或回车隔开。例：A=[123;456;789]或A=[1,2,34,5,67,8,9]2）利用函数建立很多函数都会返回一个或几个矩阵作为计算结果。以下函数则专门用于建立特殊矩阵（向量）。llinspace：建立等差数列，x=linspace(startValue,endValue,nelements)，若省略nelements，则为100。例：>>v=linspace(0,9,4)v=0369llogspace：建立等比数列。y=logspace(a,b,n)：在10^a与10^b之间产生n个数，若省略n，则为50y=logspace(a,pi,n)：在10^a与pi之间产生n个数，若省略n，则为50。例：>>y=logspace(0,2,5)y=1.00003.162310.000031.6228100.0000特殊矩阵的建立（参见p83）：l空矩阵：[]，产生行数、列数均为0的矩阵。l对角线矩阵：diag(V)用向量V建立一个对角线矩阵，diag(A)从矩阵A中提取主对角线向量。diag(A,k)提取主对角线上方第k条对角线。l单位矩阵：eye(n)l幺矩阵：ones(m),ones(m,n)产生元素值均为1的矩阵l零矩阵：zeros(m),zeros(m,n)产生元素值均为0的矩阵l随机矩阵：元素值为随机值。rand(n)、rand(m.n)：分别建立n*n、m*n阶矩阵，矩阵元素为0~1之间均匀分布的数；randn(n)、randn(m,n)：0~1之间正态分布。3）通过mat文件保存和建立矩阵savefilenamevar1var2...：把当前工作空间(workpalce)中的指定变量保存到文件filename（扩展名为.mat）中，如果不指定变量在，则保存所有变量。loadfilenamevar1var2...：把文件filename.mat中的变量载入到当前工作空间中。如果不指定变量在，则载入所有变量。也通过workplacebrowser的工具栏按钮可以完成以上操作。4)通过workplacebrowser编辑矩阵双击workplacebrowser中的变量，或选中变量后单击工具栏的open按钮，可以打开ArrayEditor编辑矩阵。\n54一、矩阵的基本操作1、矩阵元素的操作：如A(3,2)=200，如果给出的行数或列数大于原来矩阵的]范围，则会自动扩展原来的矩阵，并将扩展后未赋值的元素置0。2、冒号操作：冒号是一个重要的运算符，利用冒号可以产生向量、拆分矩阵。产生向量的一般格式是：e1:e2:e3。其中e1为初始值；e2为步长，若省略则为1；e3为终止值。例：>>t=0:2:10t=0246810拆分矩阵：A(i,:)提取A矩阵的第i行元素。A(:,j)提取A矩阵的第j列元素A(i:i+m,k:k+m)：提取A矩阵的第i—i+m行元素、第k—k+m列元素。例：A=[12345;678910;1112131415;1617181920];A(3,:)ans=1112131415A(2:3,4:5)ans=9101415A(2:3,1:2:5)ans=6810111315>>ans(:)ans=61181310153、矩阵的合并大矩阵可由方括号中的小矩阵建立起来。例：A=[123;456;789];C=[Aeye(size(A))]C=123100456010\n547890014、删除矩阵元素利用空矩阵，可将矩阵的某些行和列删除。如A(:,[2,4])=[]将删除A的2列和第4列元素。例：A=[12345;678910;1112131415];A(:,[2,4])=[]%删除A的2列和第4列元素A=13568101113155、矩阵的加减法运算规则：对应元素相加、减。A+B，A-B6、数乘：k*A7、矩阵的乘法运算规则：按线性代数中矩阵乘法运算进行，即放在前面的矩阵的各行元素，分别与放在后面的矩阵的各列元素对应相乘并相加。即例：A=[123;456];B=[12;30;74];A*Bans=281461328、数组乘法：将矩阵的对应元素相乘。x.*y例：x=[12;34];y=[56;78];z=x.*yz=51221329、矩阵的除法有两种除法运算：\（左除）和/（右除）。A\B相当于inv(A)*B，B/A相当于B*inv(A)。对于标量来说，/就是通常所用的除法，且A/B和B\A是等价的。矩阵的除法在讲线性方程组求解再详细讲解。10、数组除法：将矩阵的对应元素相除。x.\y和y./x是等价的。11、矩阵的乘方：^例：A=[12;34];A^2ans=710152212、数组的乘方：求矩阵的每个元素的幂，.^例：A=[12;34];A.^2ans=14\n5491612、转置把A的所有列换成相应的列，得到的矩阵称为A的转置矩阵。运算符是'。（对称矩阵A'=A；反对称矩阵A'=-A）例：A=[12;34];A'ans=13241.1表达式算术表达式：+,-,*,/,\,^（乘方）关系表达式：<,<=,>,>=,==,~=，若结果成立，关系表达式结果为1，否则为0。逻辑表达式：&（与），|（或），~（非）。非零元素作为真值，结果为真时表示为1，否则为0。运算优先级：算术运算优先级最高，逻辑运算最低。1.2数学函数1.2.1常用数学函数参见p25-261.2.2矩阵的超越函数一、矩阵指数Y=expm(X)，二、矩阵对数Y=logm(X)通常有：logm(expm(X))=X=expm(logm(X))三、矩阵平方根Y=sqrtm(X)习题：教材p31-32：2-(1)(2)(3)、3、4题\n54Matlab程序设计1.1m文件matlab有两种工作方式，一种是交互式的命令行工作方式，一种是m文件的程序工作方式。m文件是后缀名为.m的文件，m文件有两类：命令文件和函数文件，命令文件没有输入参数，也不返回输出参数。函数文件可以有输入参数，也可返回输出参数。函数文件和命令文件的执行如同普通的matlab命令。当输入文件名时(如果有自变量就一起附带参数，就执行文件中的语句。一、m文件的建立与编辑所有的M文件都是普通的ASCII文件，可用任何文本编辑器建立和编辑，但一般用matlab自带的编辑器，在其中可以设置程序断点，进行调试。m文件的建立：File->New->m-filem文件的建立：File->Open…二、命令文件命令文件又称脚本文件（scriptfile），它是将多条需要运行的命令编辑保存到一个文件中，然后在命令窗口中输入文件名，就可运行文件中包含的命令序列。matlab将首先在当前工作目录下寻找此文件，如果它不在当前目录下，那么在该路径下的所有目录中搜索。搜索路径可以通过File->setpath…菜单指定。例：p341.2数据的输入输出matlab提供了一些输入输出函数，便于用户和计算机之间进行交互。1.ipput函数A=input(提示信息,选项)，选项可省略，为's'时允许用户输入字符串。如：>>A=input('pleaseinputA:')pleaseinputA:[36;15]A=\n5436152.pause函数pause(延迟秒数)，省略秒时，需要按任意键后才继续执行3.disp函数disp(输出项)，输出项可以是字符串或矩阵。1.1程序结构一、选择结构1.if语句if语句根据逻辑表达式的值，决定是否执行相应的语句。其格式为：ifexpression1commandsevaluatedifexpression1isTrueelseifexpression2commandsevaluatedifexpression2isTrueelseifexpression3commandsevaluatedifexpression3isTrueelsecommandsevaluatedifnootherexpressionisTrueend如果在表达式中的所有元素为真(非零)，那么就执行if和end语言之间的{commands}。在表达式包含有几个逻辑子表达式时，即使前一个子表达式决定了表达式的最后逻辑状态，仍要计算所有的子表达式。2.switch语句switch语句根据表达式取值的不同，分别执行不同的语句。其格式为：switchexpression(scalarorstring)casevalue1statements%Executesifexpressionisvalue1casevalue2statements%Executesifexpressionisvalue2…otherwisestatements%Executesifexpressiondoesnotmatchanycase二、循环1.for语句For循环允许一组命令以固定的和预定的次数重复。For循环的一般形式是：forx=array{commands}end\n54在for和end语句之间的{commands}按数组中的每一列执行一次。在每一次迭代中，x被指定为数组的下一列，即在第n次循环中，x=array(:,n)。2.while语句与For循环以固定次数求一组命令的值相反，While循环以不定的次数求一组语句的值。While循环的一般形式是：whileexpression{commands}end只要在表达式里的所有元素为真，就执行while和end语句之间的{commands}。通常，表达式的求值给出一个标量值，但数组值也同样有效。在数组情况下，所得到数组的所有元素必须都为真。3.continue语句：跳过后面的语句，进入下一次迭代。4.break语句：跳出循环1.1函数文件主函数与子函数：除内联函数外的所有函数都需要在m文件中定义。每个M文件中可以定义一个或多个函数，最先出现在该m文件中的为主函数，其余的为子函数。主函数的范围比子函数要广，可以在M文件外部调用（在命令窗口或者是其他的M文件中），而子函数则只在主函数和该M文件的其他子函数中可见。私有函数：私有函数是主函数文件的一种。私有函数放在以private命名的子目录下，他们只是对其父目录中的函数是可见的。因此他们可以采用与其他目录下函数相同的名字，这在当你想创建自己特定的函数的新版本而想在另外目录保存原来版本的函数的时候很有用。全局变量和局部变量：函数文件与命令文件在通信方面是不同的，函数内部的变量是局部变量，不出现在MATLAB工作空间。但在函数文件中可以用global定义全局变量，全局变量的作用域是整个matlab工作空间。例：globala，则a为全局变量。全局变量应尽量少使用。一、函数文件格式function输出形参表=函数名(输入参数表)%注释说明部分函数体二、函数文件的规则和属性1.函数名和文件名必须相同。例如，函数fliplr存储在名为fliplr.m文件中。2.在函数M文件中，到第一个非注释行为止的注释行是帮助文本。当需要帮助时，返回该文本。例如，»helpfliplr返回fliplr的注释部分。第一行帮助行，名为H1行，是由lookfor命令搜索的行。3.函数可以有零个或更多个输入参量。函数可以有零个或更多个输出参量。4.从函数M文件内可以调用命令文件。在这种情况下，命令文件查看函数工作空间，不查看MATLAB工作空间。5.函数可以递归调用。即M文件函数能调用它们本身。例：p42，求n!。在编制要递归调用的函数时，必须确保会终止，否则MATLAB会陷入死循环。6.当函数M文件到达M文件终点，或者碰到返回命令return，就结束执行和返回。return命令提供了一种结束一个函数的简单方法，而不必到达文件的终点。\n547.MATLAB函数error在命令窗口显示一个字符串，放弃函数执行，把控制权返回给键盘。这个函数对提示函数使用不当很有用，如在以下文件片段中：iflength(val)>1error('VALmustbeascalar.')end这里，如果变量val不是一个标量，error显示消息字符串，把控制权返回给命令窗口和键盘。一、内联函数内联函数用于定义简单的函数，他可以直接在命令窗口中定义。内联函数的定义为：g=inline('expr','arg1','arg2',...)。例：g=inline('sin(alpha*x)','x','alpha');g(0.5,pi)ans=1二、函数调用格式：[输出参数表]=函数名[输入参数表]函数调用的注意事项：1.函数调用时，各实参出现的顺序要与函数定义时形参的顺序一致。函数可以按少于函数文件中所规定的输入和输出变量进行调用，但不能用多于函数M文件中所规定的输入和输出变量数目。2.当函数有一个以上输出变量时，输出变量包含在括号内。例如，[V,D]=eig(A)。不要把这个句法与等号右边的[V,D]相混淆。右边的[V,D]是由数组V和D所组成。3.利用系统变量nargin和nargout，可以得到调用时的输入参数和输出参数的数目，从而决定函数如何处理。nargin为输入参量个数；nargout为输出参量个数。例如，考虑MATLAB函数linspace：functiony=linspace(d1,d2,n)ifnargin==2n=100;endy=[d1+(0:n-2)*(d2-d1)/(n-1)d2];如果用户只用两个输入参量调用linspace，例如linspace(0,10)，linspace产生100个数据点。相反，如果输入参量的个数是3，例如linspace(0,10,50)，第三个参量决定数据点的个数。3.函数句柄有些函数的输入参数也是函数（如fplot），这时可以通过函数句柄来传入。函数句柄的定义方法是：在函数名前加一个@。例：fhandle=@sin;fplot(fhandle,[0,2*pi])可以使用格式：[y1,y2,...]=feval(fhandle,x1,...,xn)来调用函数句柄fhandle所指向的函数。例：x=feval(@sin,pi/2)x=1习题：教材p49-51：2、3、4、5、6题\n54图形与声音1.1二维图形一、plot函数plot(x,y,'属性值')：plot(x,y1,'属性值'，x,y2,'属性值'…)h=plot(…)以x为横坐标，y为纵坐标，按照坐标(xj,yj)的有序排列绘制曲线。属性值用于定义图形属性，如线型、颜色等（见p61），属性值可省略。例：>>x=-pi:.1:pi;y1=sin(x);y2=cos(x);plot(x,y1,'b-',x,y2,'r')二、fplot函数fplot绘制函数曲线，它可以自适应地对函数采样。格式：fplot(函数,变量取值范围,相对误差,选项)绘制函数的可以是标准函数、表达式，或用户自定义的函数。函数形式必须为y=f(x)。变量取值范围为[xminxmax]，或[xminxmaxyminymax]，默认的相对误差为2e-3。例：fplot('[sin(x),cos(x)]',[0,2*pi],1e-3,'*')fplot(@sin,[0,2*pi])%使用函数句柄fplot(inline('sin(x)'),[0,2*pi])%使用内联函数三、特殊坐标图形双对数坐标系：loglog(x,y)。两个坐标轴均用以10底的对数刻度标定。单对数坐标系：semilogx(x,y),x轴用以10底的对数刻度标定。semilogy(x,y)，y轴用以10底的对数刻度标定极坐标：polar(theta,r)，向量theta的元素代表弧度参数，向量r代表从极点开始的长度。四、其他图形函数stairs(x,y)：阶梯图形bar(x,y)：条状图形fill(x,y,'c')：填充图形。绘制并填充二维多边图形，x和y为多边形顶点向量，c规定填充颜色。errorbar(x,y,e)：误差条形图。绘制向量y对x的误差条形图，误差条对称地分布在y的上方和下方，长度为e。1.2三维图形一、plot3函数plot3(x,y,z)：用(xi,yi,\n54zi)所定义的点绘制三维曲线。向量x、y和z必须为等长度的。该命令与plot类似。例：绘制三维螺旋线t=0:pi/50:10*pi;plot3(sin(t),cos(t),t)一、mesh函数用法：h=mesh(X,Y,Z,C)功能：生成由X，Y和Z指定的网线面。若X与Y均为向量，length(X)=n，length(Y)=m，而[m，n]=size(Z)，则空间中的点(X(j),Y(I),Z(i,j))为所画曲面网线的交点，分别地，X对应于Z的列，Y对应于Z的行。C是一个矩阵，用于指定三维网格点的颜色。若省略C，则颜色由Z指定。h为返回的图形对象句柄，可省略。例：>>[x,y]=meshgrid(0:0.15:2*pi,0:0.15:2*pi);z=sin(y).*cos(x);mesh(x,y,z)(a)mesh图(b)surf图(c)waterfall图二、surf函数用法：h=surf(X,Y,Z,C)功能：绘制三维曲面图。将网格线之间的补面(patch)用颜色填充。例：接上个例子，将mesh(x,y,z)改为surf(x,y,z)。三、waterfall函数用法：h=waterfall(X,Y,Z)功能：绘制瀑布图。例：接上个例子，将mesh(x,y,z)改为waterfall(x,y,z)。1.2图形窗口控制figure：创建一个新图形窗口figure(h)：激活图形窗hclose(h)：关闭图形窗hcloseall：关闭所有图形窗口holdon：设置图形保持状态holdoff：关闭图形保持子图：subplot(m,n,p)将图形窗口分割成m行n列，并设置p所指定的子窗口为当前窗口。子窗口按行由左至右，由上至下进行编号。\n541.1图形控制一、图形标记title：Titlesfor2-Dand3-Dplotsxlabel,ylabel,zlabel：labelsfor2-Dand3-Dplotstext(x,y,'string')：addsthestringinquotestothelocationspecifiedbythepoint(x,y).gtext('string')：Placetexton2-Dgraphusingmouselegend：Graphlegendforlinesandpatches。添加图例texlabel：ProducetheTeXformatfromcharacterstring二、坐标轴设置axis函数设置坐标轴的，其格式较多，参见p63。如axis([xminxmaxyminymax])，设置坐标轴的最大值和最小值。axisoff，关闭坐标轴。三、观察点view(az,el)设置观察点角度。标量az是方位角，即观察点在xy平面上投影得到的角度。el为仰角，表示观察点在xy平面上下方的角度。az和el都用度来衡量。以上内容都可以在Figure窗口中，通过菜单或工具栏按钮进行操作。1.2动画制作动画的一般方法为：利用循环把多个图形保存到一个列向量中，然后再播放出来。fori=1:kplottingcommandsM(i)=getframe;%拷贝一个图形窗口到一个列向量中。endmovie(M,n,fps)%播放动画，n为播放次数，缺省为1；fps为帧数/秒，缺省为12。例：播放直径不断变化的球[x,y,z]=sphere;fori=1:30surf(i*x,i*y,i*z);M(i)=getframe;endmovie(M,3)%Playthemovie3times1.3声音lsound：将向量y传送给扬声器。向量确定了最大振幅。支持大多数unix平台。格式：sound(y,Fs)，Fs为采样频率，缺省为8192Hzlwavplay：PlaysoundonPC-basedaudiooutputdevice格式：wavplay(y,Fs)，Fs为采样频率，缺省为11025Hz\n54lwavrecord：RecordsoundusingPC-basedaudioinputdevice格式：y=wavrecord(n,Fs)，Fs为采样频率，缺省为11025Hzlwavread：ReadMicrosoftWAVE(.wav)soundfile格式：[y,Fs,bits]=wavread('filename',N)，Fs为采样频率，bits为采样位数（8或16），N为读取的采样点数。lwavwrite：WriteMicrosoftWAVE(.wav)soundfile格式：wavwrite(y,Fs,N,'filename')例1：x=sin(linspace(0,10000,10000));%一个纯正弦波。sound(x);%产生声音。例2：Fs=11025;y=wavrecord(2*Fs,Fs);wavplay(y,Fs);习题：教材p81-42：1、2、3、4题\n54线性代数1.1矩阵一、矩阵的行列式函数：d=det(A)，求得方阵A的行列式d。若det(A)=0，则称A为奇异矩阵，否则称为非奇异矩阵。算法：将A作LU分解。二、矩阵的条件数条件数（2-范数）的定义为：。矩阵A的条件数是一个不小于1的实数。函数c=cond(A)%求A的2-范数的条件数，即A的最大奇异值和最小奇异值的商。说明线性方程组AX=b的条件数用来衡量关于数据中的扰动，也就是A或b对解X的灵敏度。一个差条件的方程组的条件数很大。三、矩阵的秩矩阵A中最高阶非零子式的阶数称r为A的秩(rank)。矩阵的秩代表矩阵中独立方程式个数。如果一矩阵的秩数和其矩阵的列数相等，则此矩阵为满秩矩阵，否则称为降满秩矩阵。matlab函数：k=rank(A)例：>>A=[2-1-1111-2-14-62-236-97];>>rank(A)ans=3四、矩阵的初等变换下面三种变换称为矩阵的初等行变换：1）互换两行；2）以一个非零数乘以某一行；3）把某一行的k倍加到另一行上。若将定义中的“行”换成“列”，则称之为初等列变换，初等行变换和初等列变换统称为初等变换。若矩阵A经有限次初等变换变成矩阵B，则称A与B等价，记作A<->B。矩阵A作一次初等行变换后得到的矩阵B等于以一个相应的初等阵P左乘A，即B=PA矩阵A作一次初等列变换后得到的矩阵C等于以一个相应的初等阵P右乘A，即C=AP初等变换不改变矩阵的秩。\n54对秩为r的矩阵A，必可经有限次初等行变换，使A化为阶梯形；经有限次初等变换，使A化为标准形，即，其中Er为r阶单位矩阵。matlab函数：R=rref(A)功能：用初等行变换将矩阵A变为行简化阶梯形（Reducedrowecholonform）R。其每行第一个非0元素均为1。例：接上例>>R=rref(A)R=10-1001-1000010000从而可以看出A的秩为3。一、矩阵求逆若方阵A、B满足AB=BA=E，则称A、B是可逆的，B是A的逆矩阵，记为B=A-1。一些重要结论：l可逆矩阵A的逆矩阵是唯一的；ln阶矩阵A可逆的充分必要条件是：det(A)<>0，或rank(A)=n（满秩）；否则为不可逆矩阵。l若A，B可逆，则A-1、AB、A'均可逆且(A-1)-1=A，(AB)-1=B-1A-1，(A')-1=(A-1)'。逆矩阵的求法公式法：A-1=A*/det(A)。A*称为A的伴随矩阵，它是A的每个元素换成其对应的代数余子式，然后再转置得到。初等行变换法：任何满秩矩阵都能经过一系列初等行变换化为单位矩阵E，同时用一系列同样的初等行变换作用到E上，最终I就化成了A-1。即：(P1P2…Pm)(C|E)=(E|A-1)matlab函数：B=inv(A)，求方阵A的逆矩阵B。若A为奇异阵或近似奇异阵，将给出警告信息。例；接上例，B=inv(A)，会提示Warning:Matrixissingulartoworkingprecision.1.2向量空间一、向量的线性相关性对一组n个向量v1,v2,…,vn，若存在不全为0的数k1,k2,…,kn，使成立：k1v1+k2v2+…+knVn=0，称这n个向量是线性相关的，否则是线性无关的。若向量组中r个向量线性无关，而任意r+1个向量线性相关，则称向量组的秩为r。若向量组的秩为r，则向量组中任意r个线性无关的向量都称为向量组的极大无关组。矩阵的秩等于其行向量组的秩，也等于其列向量组的秩。因此可用rank()函数求向量组的秩。以列向量构建矩阵A，把A用初等行变换化为阶梯形矩阵B，则根据B的主元可求出极大无关组。例：5.1.2的例子中，>>[R,jb]=rref(A)R=10-1001-10\n5400010000jb=124jb表明第1、2、4列向量，即A(:,jb)是一个极大无关组（主元在1、2、4列上），。一、向量空间的基与维若向量空间V中的任意一个向量，都可以由一组线性无关的向量a1,a2,…,ar线性表出，则称向量组a1,a2,…,ar为V的一个基，其向量的个数r称为V的维数。向量空间的基就是向量空间这个向量组的极大无关组，维数就是这个向量组的秩。二、向量的内积与范数1.向量的内积（点积）若a=[a1,a2,…an]，b=[b1,b2,…,bn]，则a,b的内积为：=a1b1+a2b2+…+anbn=ab’matlab函数：C=dot(A,B)。若A,B为行向量，相当于A*B’，或sum(A.*B)2.向量的范数（长度）向量的2范数：。若||a||=1，则称a为单位向量。对3维向量，上式表示空间坐标系下，起点为原点，终点为（v1,v2,v3）的有向线段的长度。Matlab计算范数：n=norm(x)，返回x的2范数n。3.向量的夹角向量a,b的夹角为：cosθ=/(||a||||b||)，若θ=0，则称a,b正交。在空间坐标系中，由于=(||a||||b||)cosθ，所以a,b的内积为a的长度与b在a上的投影长度值的乘积，物理上力所作的功就是这一概念。Matlab求向量夹角：theta=subspace(A,B)，A,B为列向量例：求力G，N的合力F的大小和方向。G=[0-100];N=[50,50];F=G+N%F=[50-50]Fh=norm(F)%F的大小，70.7017theta=subspace(F’,G’)*180/pi%F与G的夹角，45度。或theta=acos(F*G'/(norm(G)*norm(F)))*180/pi4.向量叉乘两向量的叉乘是一个过两相交向量的交点，且垂直于两向量所在平面的向量。Matlab中函数：C=cross(A,B)%若A、B为向量，则返回A与B的叉乘，即C=A×B，A、B必须是3个元素的向量；若A、B为矩阵，则返回一个3×n矩阵，其中的列是A与B对应列的叉积，A、B都是3×n矩阵。例：计算垂直于向量(1,2,3)和(4,5,6)的向量。>>a=[123];b=[456];c=cross(a,b)结果显示：c=-36-3，可得垂直于向量(1,2,3)和(4,5,6)的向量为±(-3,6,-3)。\n54一、标准正交基与正交矩阵若向量a1,a2,…,ar两两正交，则称为正交向量组。若将a1,a2,…,ar单位化，则得到标准正交向量组。向量空间V中由标准正交向量组构成的基称为标准正交基（规范正交基）。若方阵满足AA'=A'A=E，则称A为正交矩阵。A为正交矩阵的充要条件是它的n个列向量是两两正交的单位向量组。若A是正交矩阵，则：1）A的行列式等于1或-1。2）A-1=A'正交变换：若T为正交矩阵，则线性变换y=Tx的正交变换。正交变换保持向量的长度不变。Matlab中函数：Q=orth(A)，求矩阵A代表的列向量组的标准正交基。B'*B=eye(rank(A))。例：接5.1.2的例>>B=orth(A)B=-0.05840.37540.1697-0.10360.0946-0.98060.36660.8682-0.0018-0.92280.31050.09861.2线性方程组一、齐次线性方程组的解齐次线性方程组AX=0有非0解的充要条件是：A的列向量是线性相关的，即A的秩r>A=[1221;21-2-2;1-1-4-3];>>formatrat%指定有理式格式输出>>B=null(A,'r')%求解空间的有理基B=25/3%rank(A)=2，所以B有两列-2-4/31001%从而方程组的通解为：X=[2*k1+5/3*k2][-2*k1-4/3*k2][k1]\n54[k2]一、非齐次线性方程组的唯一解或特解非齐次线性方程组AX=b有解的充要条件是：列向量b可由A的列向量线性表出，即A的秩r1=增广矩阵[Ab]的秩r2。方程组有唯一解的充要条件是：b可由A的列向量线性表出，且表示法唯一，即A的秩r1=增广矩阵的秩r2=A的列数n。若r1=r2>A=[11-3;3-12;15-9];B=[140]';r=rank(A)%求秩X=A\B%求解X=1.2500-0.2500-0.0000也可用rref求解：>>C=[A,B];%构成增广矩阵R=rref(C)R=1.0000001.250001.00000-0.2500001.00000R的最后一列元素就是所求之解。2．利用矩阵的分解求方程组的解1）LU分解：LU分解又称Gauss消去分解，可把任意方阵分解为下三角矩阵的基本变换形式（行交换）和上三角矩阵的乘积。即A=LU，L为下三角阵，U为上三角阵。则方程A*X=b变成L*U*X=b，所以X=U\(L\b)，这样可以大大提高运算速度。\n54格式：[L，U]=lu(A)例：解上题>>[L,U]=lu(A)X=U\(L\B)L=0.33330.25001.00001.0000000.33331.00000U=3.0000-1.00002.000005.3333-9.666700-1.2500X=1.2500-0.2500-0.00002）Cholesky分解若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成：A=R'R，其中R为上三角阵，R'为下三角阵。则方程A*X=b变成R'RX=b，所以X=R\(R'\b)。格式：R=chol(A)3）QR分解对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形式，即：A=QR，则方程AX=b变成QRX=b，所以X=R\(Q\b)格式：[Q,R]=qr(A)这三种分解在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。3.线性方程组的迭代解法原理：将Ax=b改写成一个等价的方程组x=Bx+g，从而建立迭代式：xk+1=Bxk+g给出初始向量x0后，按公式进行迭代，就可以得到一个向量序列{xk}，若xk收敛于确定的向量x，则x就是方程组的解。雅可比迭代法：xk+1=D-1(E+F)xk+D-1b高斯-赛德尔迭代法：xk+1=(D-E)-1Fxk+(D-E)-1b线性方程组的LQ解法x=symmlq(A,b)双共轭梯度法解方程组x=bicg(A,b)稳定双共轭梯度方法解方程组x=bicgstab(A,b)复共轭梯度平方法解方程组x=cgs(A,b)共轭梯度的LSQR方法x=lsqr(A,b)广义最小残差法x=gmres(A,b)最小残差法解方程组x=minres(A,b)预处理共轭梯度方法x=pcg(A,b)准最小残差法解方程组x=qmr(A,b)例调用MATLAB数据文件west0479。>>loadwest0479>>A=west0479;%将数据取为系数矩阵A。>>b=sum(A,2);%将A的各行求和，构成一列向量。>>X=A\b;%用“\”求AX=b的解。\n54>>norm(b-A*X)/norm(b)%计算解的相对误差。ans=1.2454e-017>>[x,flag,relres,iter,resvec]=bicg(A,b)%用bicg函数求解。x=（全为0，由于太长，不显示出来）flag=1%表示在默认迭代次数（20次）内不收敛。relres=%相对残差relres=norm(b-A*x)/norm(b)=norm(b)/norm(b)=1。1iter=%表明解法不当，使得初始估计值0向量比后来所有迭代值都好。0一、求非齐次线性方程组的通解非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解。因此，步骤为：第一步：判断AX=b是否有解，若有解则进行第二步第二步：求AX=b的一个特解第三步：求AX=0的通解第四步：AX=b的通解=AX=0的通解+AX=b的一个特解。例求解方程组解：在Matlab中建立M文件如下：A=[1-23-1;3-15-3;2-12-2];b=[123]';B=[Ab];n=4;R_A=rank(A)R_B=rank(B)formatratifR_A==R_B&R_A==n%判断有唯一解X=A\belseifR_A==R_B&R_An，也就是说方程的个数多于未知数，则称为超定方程组。通常这些方程组是矛盾的，所以方程组没有精确解。在拟合实验数据的曲线时，常会遇到这个问题。关键是要找到一个向量x使它对m个方程的总误差最小。有几种方法可以求得，但是最常用的方法是最小二乘法。最小二乘解可以用除法运算符\或命令nnls和lscov来解得。MATLAB在用左除解超定方程组时都用QR因式分解。如果方程组的系数矩阵不满秩，就没有唯一的最小二乘解，将不定值设为零。然后给出一个警告信息。1.2特征值与二次型一、特征值与特征向量的基本概念设A为n阶矩阵，若存在复数λ0和非零n维列向量X使得AX=λX（1）则称数λ为A的一个特征值(eigenvalue)，X为A的属于（或对应于）特征值λ的特征向量(eigenvector)。（1）式可表达为：λX—AX=（λE—A）X=0（2）由此可见，特征值λ是使齐次线性方程组（λE—A）X=0有非零解X≠0的λ值。方程组（2）有非零解的充要条件是：|λE－A|=0（或秩（λE－A））。|λE－A|=0称为A的特征方程，|λE－A|称为A的特征多项式。A的特征值就是特征方程的根广义特征值问题：求方程解的问题。若A的特征值为λ1、λ2、λn，则有：λ1+λ2+…+λn=trace(A)λ1λ2…λn=det(A)二、特征值与特征向量的求法雅可比法是用来求实对称矩阵全部特征值与特征向量的方法，其基本思想是把对称阵A经一系列正交相似变换，化为一个对角阵，对角阵的对角元就是A的特征值，而由正交变换的矩阵的乘积可求得特征向量。QR方法是一种求一般矩阵的全部特征值的有效方法。其原理是：利用矩阵A的QR分解，得到A的一系列相似矩阵。设A=A1=Q1R1，令A2=R1Q1=Q2R2，…，Ak=Rk-1Qk-1=QkRk。由于Ak-1=Qk-1Rk-1，所以有Ak=Rk-1Qk-1=Q'k-1Ak-1Qk-1，即Ak相似于Ak-1，因此A1，A2，…，Ak均相似，它们具有相同的特征值。当k—>∞时，Ak的对角线下（不包括对角元）元素收敛于0，根据相似关系，Ak的对角元收敛于A的特征值。从而可以求出A的特征值，根据特征值由（λE—A）X=0可以解出相应的特征向量。格式E=eig(A)%求矩阵A的特征值E，以向量形式存放E。[V,D]=eig(A)%计算A的特征值对角阵D和特征向量V，使AV=VD成立。特征值在D的对角线上，对应的特征向量为V的列。V的列为单位向量。\n54[V,D]=eig(A,B)%计算广义特征值向量阵V和广义特征值阵D，满足AV=BVD。例：>>A=[1-22;-2-24;24-2];[V,D]=eig(A)V=0.33330.9339-0.12930.6667-0.3304-0.6681-0.66670.1365-0.7327D=-7.00000002.00000002.0000即：特征值-7对应特征向量(0.33330.6667-0.6667)T特征值2对应特征向量(0.9339-0.33040.1365)T和(-0.1293-0.6681-0.7327)T一、提高特征值的计算精度函数balance格式[T,B]=balance(A)%求相似变换矩阵T和平衡矩阵B，满足。B=balance(A)%求平衡矩阵B二、实对称矩阵的相似对角化相似矩阵：设A、B为n阶方阵，若存在n阶可逆矩阵P，使P-1AP=B，A相似于B，记作A~B。称P为A到B的相似变换矩阵。若A、B相似，则：A、B的特征值相同；Bm=(P-1AP)m=P-1AmP，因此若A为对角阵，则很容易通过P-1AmP求得Bm，这在很多工程问题中有重要作用。实对称矩阵（A'=A）一定可以相似对角化。[V,D]=eig(A)求得的V就是相似变换矩阵，D就是A的相似对角矩阵。V-1AV=D例：>>A=[324;202;423];[V,D]=eig(A);inv(V)*A*Vans=-1.00000.0000-0.00000.0000-1.00000.00000.0000-0.00008.0000求线性常微分方程组的非0解A=[010;001;-1243];[V,D]=eig(A)V=[V(1,:)./V(1,:);V(2,:)./V(1,:);V(3,:)./V(1,:)]\n54运行结果为V=0.10483-0.21822-0.218220.31449-0.436440.436440.94346-0.87287-0.87287D=30002000-2V=11132-2944即解为：，其中C1,C2,C3为任意常数。一、二次型及其标准型n个变数x1、x2…xn的二次齐次多项式f(x1,x2,…xn)称为二次型。一般将二次型写成如下形式：式中aij=aji,A是n阶[实]对称阵，并称这个实对称阵为二次型的矩阵。显然，二次型与实对称矩阵是一一对应的。若在上式中对一切i≠j均有aij=0，即二次型的矩阵A为对角阵时，称这样的二次型具有标准形。若存在n阶可逆矩阵P，使C'AC=B，则称方阵A合同于B，记作。作线性变换x=Py，则f=X'AX=(Py)'Apy=y'(P'AP)y。因此二次型f能否化为标准形等价于实对称矩阵A能否与对角阵合同。二次型的标准型不是唯一的，但其正系数与负系数的个数是不变的。对于实二次型f=x'Ax，若对任意x≠0，都有f>0，则称f为正定二次型，并称矩阵A是正定的，记作A>0。若对任意x≠0，都有f<0，则称f为负定二次型，并称矩阵A是负定的，记作A<0。实二次型f正定的充要条件可以用以下几种方法表述l二次型的标准型的n个系数全为正。l二次型的对称矩阵A的特征值全为正。l霍尔维茨定理：对称矩阵A的各阶顺序主子式全大于零，则其为正定；A的奇数阶主子式为负，偶数阶主子式为正，则其为负定。例：求一个正交变换X=PY，把二次型化成标准形。A=[1-22;-2-24;24-2];[V,D]=eig(A)V=0.333330.93391-0.129250.66667-0.33042-0.66812-0.666670.13654-0.73274\n54D=-700020002因此通过正交变换X=VY，化为标准型验证X=[13;21;37];Y=inv(V)*X;X'*A*XY'*D*Yans=3597ans=3597\n54习题1.分别用rank和rref函数求矩阵的秩。2.解矩阵方程：3.已知向量组a1=[1,4,2,1],a2=[-2,1,5,1],a3=[-1,2,4,1],a4=[-2,1,-1,1],a5=[2,3,0,1/3],求该向量组的秩，并求一个极大无关组。4.已知R4上的两个向量u=[4,1,2,3]'，v=[-2,3,1,4]'，求u,v的长度，以及u在v上的投影向量。5.将矩阵正交规范化。6.求方程组的解。7.求方程组的通解。8.求矩阵的特征值和特征向量。9.将矩阵A=[-102;012;220]相似对角化。10.将二次型化为标准型。11.解微分方程（提示：令x=x1,dx/dt=x2）\n54数据处理与多项式1.1基本统计处理l求最大值c=max(x),[c,I]=max(x)返回各列中的最大数，I为最大数所在的列。l求最小值min(x),[c,I]=min(x)例：>>v=[3-24-150];max(v)ans=5>>[vmin,kmin]=min(v)vmin=-2kmin=2l求中位数median(x)：若数据序列为奇数个数时，中位数就是中间那个数；个数为偶数时，为中间两个数的平均值。例：>>x1=[9–25712];median(x);ans=7>>x2=[9–256712];median(x);ans=6.5l升序排序sort(x)l求和sum(x)l求平均值mean(x)例：>>data1=2*rand(1,500)+2;mean(data1)ans=2.9940l求积prod(x)例：>>nfact=prod(1:4)%求24！nfact=24l求累计和cumsum(x)l求累计积cumprod(x)例：>>cumsum(1:5)ans=[1361015]>>cumprod(1:5)ans=12624120l求标准方差std(x)：s=std(X,flag)，flag=0（缺省值）时采用(1)式，flag=1时采用(2)式\n54>>x=randn(1,500);%标准正态分布std(x)ans=0.97641.1多项式一、多项式的表示多项式用向量形式表达，例如对多项式，用p=[10-2-5]表示。若已知多项式的全部根x1,x2…xn，即(x-x1)(x-x2)…(x-xn)=0，则可用poly函数建立多项式：p=poly(x)，x为根组成的向量，p为多项式系数组成的向量。如:>>p=poly([000])p=1000%p(x)=x3二、多项式求根x=roots(c)例：p=[10-2-5];%p(x)=x3-2x-5x=roots(p)ans=2.0946-1.0473+1.1359i-1.0473-1.1359i三、求多项式的值Y=polyval(A,x)例：>>p=[10-2-5];polyval(p,[579])ans=110324706即x=5,7,9时，多项式的值分别为110,324,706。四、多项式的四则运算多项式的加、减：直接对系数向量进行加、减（如果多项式的次数不同，需要把低次的多项式系数不足的高次项用0补足）。多项式乘法：w=conv(u,v)多项式除法：[q,r]=deconv(u,v)，q为商，r为余项，即u=conv(v,q)+r例：>>u=[1800-10];v=[2-13];w=conv(u,v)%求x4+8x3-10与2x3-x3+3的乘积w=215-524-2010-30%乘积为[q,r]=deconv(v,u)q=0.50004.25001.3750\n54r=000-11.3750-14.12501.1数据插值插值法是实用的数值方法，是函数逼近的重要方法。在生产和科学实验中，自变量x与因变量y的函数y=f(x)的关系式有时不能直接写出表达式，而只能得到函数在若干个点的函数值或导数值。当要求知道观测点之外的函数值时，需要估计函数值在该点的值。如何根据观测点的值，构造一个比较简单的函数y=φ(x)，使函数在观测点的值等于已知的数值或导数值。用简单函数y=φ(x)在点x处的值来估计未知函数y=f(x)在x点的值。寻找这样的函数φ(x)，办法是很多的。φ(x)可以是一个代数多项式，或是三角多项式，也可以是有理分式；φ(x)可以是任意光滑（任意阶导数连续）的函数或是分段函数。函数类的不同，自然地有不同的逼近效果。如拉格朗日n次插值多项式可写为：特别地，一次插值（线性插值）多项式为：或命令1interp1功能一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。各个参量之间的关系示意图为图2-14。图6 1数据点与插值点关系示意图格式yi=interp1(x,Y,xi)%返回插值向量yi，每一元素对应于参量xi，同时由向量x与Y的内插值决定。参量x指定数据Y的点。若Y为一矩阵，则按Y的每列计算。yi是阶数为length(xi)*size(Y,2)的输出矩阵。yi=interp1(Y,xi)%假定x=1:N，其中N为向量Y的长度，或者为矩阵Y的行数。yi=interp1(x,Y,xi,method)%用指定的算法计算插值：’nearest’：最近邻点插值；’linear’：线性插值（缺省方式）；’spline’：三次样条函数插值。’pchip’：分段三次Hermite插值。’cubic’：与’pchip’操作相同；yi=interp1(x,Y,xi,method,'extrap')%对于超出x范围的xi中的分量将执行特殊的外插值法extrap。yi=interp1(x,Y,xi,method,extrapval)%确定超出x范围的xi中的分量的外插值extrapval，其值通常取NaN或0。例：已知1900年至1990年中每隔10年的美国人口数（万）t=1900:10:1990;\n54p=[759991971057112320131661506917932203212265024963];interp1(t,p,1975)%求1975年的人口数ans=21486x=1900:1:2000;%绘制分段线性插值曲线y=interp1(t,p,x,'spline');plot(t,p,'o',x,y)图6 2一维数据插值命令2interp2功能二维数据内插值格式ZI=interp2(X,Y,Z,XI,YI)%返回矩阵ZI，其元素包含对应于参量XI与YI（可以是向量、或同型矩阵）的元素，即Zi(i,j)←[Xi(i,j),yi(i,j)]。用户可以输入行向量Xi和列向量Yi，此时输出向量Zi与矩阵meshgrid(xi,yi)是同型的。同时取决于由输入矩阵X、Y与Z确定的二维函数Z=f(X,Y)。参量X与Y必须是单调的，且相同的划分格式，就像由命令meshgrid生成的一样。若Xi与Yi中有在X与Y范围之外的点，则相应地返回nan（NotaNumber）。ZI=interp2(Z,XI,YI)%缺省地，X=1:n、Y=1:m，其中[m,n]=size(Z)。再按第一种情形进行计算。ZI=interp2(Z,n)%作n次递归计算，在Z的每两个元素之间插入它们的二维插值，这样，Z的阶数将不断增加。interp2(Z)等价于interp2(z,1)。ZI=interp2(X,Y,Z,XI,YI,method)%用指定的算法method计算二维插值：’linear’：双线性插值算法（缺省算法）；’nearest’：最临近插值；’spline’：三次样条插值；’cubic’：双三次插值。例：平板温度length=1:5;%indexforlengthofplate(i.e.,thex-dimension)with=1:3;%indexforwithofplate(i,e,,they-dimension)temps=[8281808284;7963616581;8484828586];%求长2.5、宽1.2处的温度interp2(length,with,temps,2.5,1.2)ans=76.8%绘制分段线性插值曲面[l,w]=meshgrid(1:0.2:5,1:0.2:3);%或l=[1:0.2:5];w=[1:0.2:3]'choosehigherresolutiont=interp2(length,with,temps,l,w,'cubic');%linearinterpolationmesh(l,w,t)xlabel('lengthofPlate'),ylabel('withofPlate')zlabel('DegreesCelsius')\n54图6 3二维数据内插值命令3interp3功能三维数据插值（查表）格式VI=interp3(X,Y,Z,V,XI,YI,ZI)%找出由参量X,Y,Z决定的三元函数V=V(X,Y,Z)在点（XI,YI,ZI）的值。参量XI,YI,ZI是同型阵列或向量。若向量参量XI,YI,ZI是不同长度，不同方向（行或列）的向量，这时输出参量VI与Y1,Y2,Y3为同型矩阵。其中Y1,Y2,Y3为用命令meshgrid(XI,YI,ZI)生成的同型阵列。若插值点(XI,YI,ZI)中有位于点(X,Y,Z)之外的点，则相应地返回特殊变量值NaN。VI=interp3(V,XI,YI,ZI)%缺省地，X=1:N，Y=1:M，Z=1:P，其中，[M,N,P]=size(V)，再按上面的情形计算。VI=interp3(V,n)%作n次递归计算，在V的每两个元素之间插入它们的三维插值。这样，V的阶数将不断增加。interp3(V)等价于interp3(V,1)。VI=interp3(…,method)%用指定的算法method作插值计算：‘linear’：线性插值（缺省算法）；‘cubic’：三次插值；‘spline’：三次样条插值；‘nearest’：最邻近插值。说明在所有的算法中，都要求X,Y,Z是单调且有相同的格点形式。当X,Y,Z是等距且单调时，用算法’*linear’，’*cubic’，’*nearest’，可得到快速插值。>>[x,y,z,v]=flow(20);[xx,yy,zz]=meshgrid(.1:.25:10,-3:.25:3,-3:.25:3);vv=interp3(x,y,z,v,xx,yy,zz);slice(xx,yy,zz,vv,[69.5],[12],[-2.2]);shadinginterp;colormapcool插值图形如图。图6 4三维插值图命令4spline功能三次样条数据插值格式yy=spline(x,y,xx)%该命令用三次样条插值计算出由向量x与y确定的一元函数y=f(x)\n54在点xx处的值。若参量y是一矩阵，则以y的每一列和x配对，再分别计算由它们确定的函数在点xx处的值。则yy是一阶数为length(xx)*size(y,2)的矩阵。pp=spline(x,y)%返回由向量x与y确定的分段样条多项式的系数矩阵pp，它可用于命令ppval、unmkpp的计算。例：对离散地分布在y=exp(x)sin(x)函数曲线上的数据点进行样条插值计算：>>x=[024581212.817.219.920];y=exp(x).*sin(x);xx=0:.25:20;yy=spline(x,y,xx);plot(x,y,'o',xx,yy)插值图形结果如图。图6 5三次样条插值1.1曲线拟合从几何意义上将，拟合是给定了空间中的一些点，找到一个已知形式未知参数的连续曲面来最大限度地逼近这些点；而插值是找到一个(或几个分片光滑的)连续曲面来穿过这些点。通常采用最小二乘拟合，以最小二乘多项式拟合为例说明。设有变量x和y的一组数据(xi,yi)，i=1,2…m。最小二乘多项式拟合就是要求取多项式P(x)=a0+a1x+…+anxn的系数a0,a1…an(nbasicfitting菜单，打开basicfitting对话框，在其中可以进行曲线拟合。图6 7basicfitting对话框1.1离散傅立叶变换傅立叶变换用于频域分析，即以频率为横坐标，研究信号的幅度、相位在各频率上的分布情况。傅立叶变换的定义为：F(ω)为复数，可表示为：幅度频谱：|F(ω)|与ω的关系图。相位频谱：φ(ω)与ω的关系图。为利用计算机计算傅立叶变换，需要对无限长连续信号，进行时域采样、时域截断、频域采样，周期延拓，从而得到离散傅立叶变换DFT：忽略由频谱混叠和泄漏引起的误差，离散傅立叶变换与傅立叶变换具有下列关系：\n54式中，f0=1/NTs，N为采样点数，Ts为采样周期。DFT方法计算量太大，限制了应用。1965年美国学者提出了一种快速计算DFT的算法,称为快速傅立叶变换FFT。matlab中的函数为：Y=fft(X,n)，n表示只计算x的前n点，省略时计算全部点。例：信号x(t)=sin(2πf1)+2sin(2πf2)+N。式中，f1=50，f2=120，N为正态分布的噪声信号。t=0:0.001:0.6;%600点x=sin(2*pi*50*t)+2*sin(2*pi*120*t)+2*randn(size(t));subplot(2,1,1)plot(x(1:256))%绘制信号的时域波形y=fft(x,512);%快速傅立叶变换f=(0:256)/(0.001*512);%频率轴，只绘制前面257点subplot(2,1,2)plot(f,abs(y(1:257)))%绘制信号的幅度频谱信号的时域波形和幅度频谱如下。从时域波形只能看出信号的振幅范围，很难看出信号的频率特征；而从频谱图可以很容易看出信号中主要由频率分别为50与120的简谐信号构成，且后者的幅度是前者的两倍。图6 8(a)信号的时域波形(b)信号的幅度频谱\n54微积分1.1数值微分函数f(x)在x=a的微分可用差分近似表示。前向差分：后向差分：中间差分：在MATLAB用diff函数来计算二相邻点的差值，语法为d=diff(y)。其中y为一系列的数据向量，d为差分值输出。d=diff(y,n)计算n阶差分。函数y=f(x)的导数则为dy=diff(y)./diff(x)。>>x=[1491625];>>dx=diff(x)dx=3579>>dx2=diff(x,2)dx2=222现假设有一组流率与压降之关系的实验数据，如下表:x(流率,L/s)010.8016.0322.9132.5636.7639.8843.68y(压降,KPa)00.2990.5761.0361.7812.4322.8463.304现利用diff来进行後向差分,藉以求得压降对流率之变化率,并与中间差分之微分近似值结果做一比较.>>x=[010.8016.0322.9132.5636.7639.8843.68];>>y=[00.2990.5761.0361.7812.4322.8463.304];>>n=length(x);>>d1=diff(y)./diff(x);%前、后向差分>>d2=0.5*(d1(1:n-2)+d1(2:n-1));%中间差分>>plot(x(2:n),d1,'o--',x(1:n-1),d1,'*-.',x(2:n-1),d2,'s-')>>xlabel('x')>>ylabel('微分近似值')>>legend('后向差分','前向差分','中间差分')1.2数值积分函数f(x)在区间[a,b]上的定积分为：S=\n54定积分的数值解法的基本思想是：将积分区间[a,b]分成n个子区间，再在每个子区间上求定积分的近似值。即。图7 1数值积分梯形求积：在每个子区间上用线性插值代替f(x)，从而得到S的数值解为：，h=(b-a)/n，称为步长。梯形求积公式的误差为：辛普生求积(Simplesonquadrature)：在每两个相邻的子区间上用二次插值代替f(x)，从而得到S的数值解为：，其误差为：实际应用时，往往采用逐步缩小步长的办法，直到两次的计算结果相差不大时为止，称为变步长积分法。缩小步长最好每次减半，这样前次计算过的函数值不必重复计算。一、一元函数的数值积分函数1quad、quadl、quad8功能求函数的数值积分。quad采用变步长Simpleson积分法，quad8采用高阶方法，新版本中quad8用quadl取代，quadl采用变步长Gauss/Lobatto方法。格式q=quad(fun,a,b)%近似地从a到b计算函数fun的数值积分，误差为10-6。若给fun输入向量x，应返回向量y，即fun是一单值函数。q=quad(fun,a,b,tol)%用指定的绝对误差tol代替缺省误差。tol越大，函数计算的次数越少，速度越快，但结果精度变小。q=quad(fun,a,b,tol,trace,p1,p2,…)%将可选参数p1,p2,…等传递给函数fun(x,p1,p2,…)，再作数值积分。若tol=[]或trace=[]，则用缺省值进行计算。[q,n]=quad(fun,a,b,…)%同时返回函数计算的次数n…=quadl(fun,a,b,…)%用高精度进行计算，效率可能比quad更好。例：求quad('sin(x)',0,pi)ans=2.0000函数2trapz功能梯形法数值积分\n54格式T=trapz(Y)%用等距梯形法近似计算Y的积分。若Y是一向量，则trapz(Y)为Y的积分；若Y是一矩阵，则trapz(Y)为Y的每一列的积分；若Y是一多维阵列，则trapz(Y)沿着Y的第一个非单元集的方向进行计算。T=trapz(X,Y)%用梯形法计算Y在X点上的积分。若X为一列向量，Y为矩阵，且size(Y,1)=length(X)，则trapz(X,Y)通过Y的第一个非单元集方向进行计算。例：求>>X=0:pi/100:pi;Y=sin(X);Z=trapz(X,Y)Z=1.9998一、重积分的数值计算函数1dblquad功能矩形区域上的二重积分的数值计算格式q=dblquad(fun,xmin,xmax,ymin,ymax)%调用函数quad在区域[xmin,xmax,ymin,ymax]上计算二元函数z=f(x,y)的二重积分。输入向量x，标量y，则f(x,y)必须返回一用于积分的向量。q=dblquad(fun,xmin,xmax,ymin,ymax,tol)%用指定的精度tol代替缺省精度10-6，再进行计算。q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method)%用指定的算法method代替缺省算法quad。method的取值有@quadl或用户指定的、与命令quad与quadl有相同调用次序的函数句柄。q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method,p1,p2,…)%将可选参数p1,p2,..等传递给函数fun(x,y,p1,p2,…)。若tol=[]，method=[]，则使用缺省精度和算法quad。例>>fun=inline(’y./sin(x)+x.*exp(y)’);>>Q=dblquad(fun,1,3,5,7)Q=3.8319e+003函数2quad2dggen（在NIT工具箱中）功能任意区域上二元函数的数值积分。格式q=quad2dggen(fun,xlower,xupper,ymin,ymax)%在由[xlower,xupper,ymin,ymax]指定的区域上计算二元函数z=f(x,y)的二重积分。例计算单位圆域上的积分：先把二重积分转化为二次积分的形式：f=inline(’exp(-x.^2/2).*sin(x.^2+y)’,’x’,’y’);xlower=inline(’-sqrt(1-y.^2)’,’y’);xupper=inline(’sqrt(1-y.^2)’,’y’);Q=quad2dggen(fun,xlower,xupper,-1,1,1e-4)计算结果为：Q=0.5368603818函数3：triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol)：长方体区域的三重积分函数4：quadndg('Fun',xlow,xhigh,tol)，在NIT工具箱中，n-Dhyper-rectangular区域的n重积分。\n541.1常微分方程的解微分方程是研究函数变化规律的重要工具，有着非常广泛的应用。若微分方程的自变量个数为1则称为常微分方程（ODE，OrdinaryDifferentialEquation）；若自变量的个数大于1，则称为偏微分方程。一般的，n阶常微分方程具有形式：若式子左端为x及的一次有理整式，则称为n阶线性微分方程。一般的，n阶线性常微分方程具有形式：高阶线性微分方程一般可化为一阶线性微分方程组，方法为：令xn=dxn-1/dt,xn-1=dxn-2/dt,…,x2=dx/dt,x1=x，则高阶的方程（组）可写成一阶微分方程组：X'=AX+f(t)。若上式中f(t)=0，则称为齐线性方程组。n阶齐线性方程组X'=AX一定存在n个线性无关的解：x1(t),x2(t),…xn(t)，且其通解可表为：x=C1x1(t)+C2x2(t)+…+Cnxn(t)，其中Ci为任意常数。x1(t),x2(t),…xn(t)称为方程组的基础解系，其构成的矩阵称为基解矩阵。非齐线性方程组的通解=对应齐线性方程组的通解+非齐线性方程组的一个解。若式中A为常数，则称为常系数线性微分方程组。为了确定微分方程的特定解，需要给出这个解必须满足的条件，即定解条件。如果定解条件是初始条件，就称为初值问题（IVP,InitialValueProblem）。一、常系数齐线性微分方程组的代数解法定理：矩阵eAt是常系数齐线性微分方程组X'=AX的基解矩阵。，如果A是一个对角矩阵，则，式中a1,a2…为A的对角线元素。若A是实的，则eAt也是实的定理：如果A具有n个线性无关的特征向量v1,v2,…vn，它们对应的特征值分别为λ1,λ2,…λn（不必各不相同），则矩阵是齐线性微分方程组X'=AX的一个基解矩阵。关于这一结论，在第5章举例说明过。如果A是实的，则通过eAt=Φ(t)Φ-1(0)，可以得到一个实的基解矩阵。如果A的特征向量并非是线性无关的，则不能通过相似变换V-1AV化成对角矩阵，而只能化成准对角矩阵或约当标准型（JordanCanonicalForm）矩阵，利用约当标准型矩阵可计算基解矩阵。\n54约当标准型矩阵，其中称为约当块，为nj阶矩阵，且。λ1，λ2，…λk是A的特征值，nj是λj的重特征值数。从而基解矩阵，而，其中matlab符号工具箱中的函数jordan(A)可以求出非奇异的矩阵V，使得inv(V)*A*V成为约当标准型，函数jordan的形式为[V,J]=jordan(A)，V为相似变换矩阵，J为约当标准型矩阵。例：>>A=[21;-14];[V,J]=jordan(A)%findtheJordanformandeigenvectorsV=-11-10J=3103iV=inv(V)iV=0-11-1从而基解矩也可利用符号工具箱求得：>>symstV*[1t;01]*exp(3*t)*inv(V)ans=[exp(3*t)*(-t+1),exp(3*t)-exp(3*t)*(-t+1)][-exp(3*t)*t,exp(3*t)+exp(3*t)*t]\n54一、利用拉普拉斯变换解齐线性微分方程组（略）二、常微分方程的数值解一个一阶常微分方程式(ordinarydifferentialequation,ODE)可以下式表示其中t为独立变数而y是t的函数。数值解法就是要求方程的解在点t0>globalMU>>MU=7;>>Y0=[1;0]>>[t,x]=ode45(‘verderpol’,0,40,Y0);>>x1=x(:,1);x2=x(:,2);>>plot(t,x1,t,x2)图7 3VerderPol微分方程图1.1非线性方程(组)求解一、非线性方程的解非线性方程的标准形式为f(x)=0求方程实根的最简单的近似方法是二分法，但其求解时间长，一般不单独采用。求方程根的主要方法是迭代法。迭代法的原理如下：将方程f(x)=0化为一个等价（同根）的方程：x=ψ(x)。给定一个初值x0，代入右端可算得一个x1=ψ(x0)，再将x1代入右端，又可算得x2=ψ(x1)，如此继续下去，得到一个序列{xk}，其中xk+1=ψ(xk)，{xk}称为迭代序列，ψ(x)称为迭代函数。如果迭代序列是收敛的，且收敛于x*，则x*就是方程的根。实际计算时，当两次迭代的值相差很小时，就可取xk+1作为方程的根。如何构造迭代函数，才能保证迭代序列收敛呢？可以将f(x)在xk处作泰勒展开，并取近似值，得f(x)=f(xk)+f'(xk)(x-xk)，从而得到近似的方程f(xk)+f'(xk)(x-xk)=0。所以函数fzero说明该函数采用数值解求方程f(x)=0的根。\n54格式x=fzero(fun,x0)%用fun定义表达式f(x)，x0为初始解。[x,fval]=fzero(fun,x0)%fval=f(x)例：求x3-2x-5=0的根>>fun='x^3-2*x-5';z=fzero(fun,2)%初始估计值为2z=2.0946[x,fval]=fzero(fun,2)x=2.0946fval=-8.8818e-016一、非线性方程组的解非线性方程组的标准形式为：F(X)=0，其中X为向量，F(X)为函数向量。求解非线性方程组常用办法是迭代法，把方程组化为如下的等价形式：li(x1,x2,…xn)=ψ(x1,x2,,…,xn)。其中左端是线性式，建立迭代式：取初值()后，按上式迭代，右端成为已知量，因而每步迭代需要求解一个线代数方程组，解这个线代数方程组若用迭代法，就形成一个双重迭代。函数fsolve格式x=fsolve(fun,x0)%用fun定义向量函数，定义方式为：先定义方程函数functionF=myfun(x)。F=[表达式1；表达式2；…表达式m]%保存为myfun.m，并用下面方式调用：x=fsolve(@myfun,x0)，x0为初始估计值。x=fsolve(fun,x0,options)[x,fval]=fsolve(…)%fval=F(x)，即函数值向量例:求下列系统的根解：化为标准形式，设初值点为x0=[-5-5]。先建立方程函数文件，并保存为myfun.m：functionF=myfun(x)F=[2*x(1)-x(2)-exp(-x(1));-x(1)+2*x(2)-exp(-x(2))];然后调用优化程序x0=[-5;-5];%初始点[x,fval]=fsolve(@myfun,x0)x=0.5671\n540.5671fval=1.0e-006*-0.4059-0.40591.1函数优化最优化是求最优解，也就是在某个区间内有条件约束或者无条件约束地找到函数的最大值或者最小值。MATLAB使用数字方法求函数的最小值。使用迭代算法，也就是有些步骤要重复许多次。现在，假设要求函数f在某个区间内的最小值xmin。迭代方法需要一个初始估计值x0。从x0开始找到一个更接近xmin的新值x1，这个值的好坏取决于使用的数学方法。直到找到有足够精度的近似值xi才停止迭代，也就是绝对值|xmin－xi|足够小。这里提到了标准MATLAB系统的两个最优化命令，fmin命令可以求单变量函数的最小值；fmins命令可以求多变量函数的最小值，同时它还要求有一个初始向量。在新版本中，fmin和fmins分别被fminbnd和fminsearch取代。Matlab中没有求函数f的最大值的命令，但可以通过求其相反函数h=－f的最小值间接求得。函数格式：x=fmin(fcn,x1,x2,options)求函数在区间(x1，x2)内取最小值时的x值，采用黄金分割法和二次插值。fcn是目标函数名。如果没有局部最小值，则返回区间内的最小x值。向量options为控制参数，如options(1)=1，显示中间结果；options(2)表示得到的结果x的精度，缺省为10－4，options可省略。x=fmins(fcn,x0,options)求函数fcn的最小值。由用户自己给出一个初始估计向量x0，优化参数options可省略。例：在区间[0，2]内求函数f(x)=x3-2x-5的最小值。f=inline('x.^3-2*x-5');x=fmin(f,0,2)x=0.8165f(x)ans=-6.0887例：求香蕉函数的最小值。>>banana=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2');x=fmins(banana,[-1.2,1])x=11banana([1,1])ans=0此外，MATLAB提供了最优化工具箱。\n54符号计算符号数学工具箱是操作和解决符号表达式的符号数学工具箱(函数)集合，有复合、简化、微分、积分以及求解代数方程和微分方程的工具。另外还有一些用于线性代数的工具，求解逆、行列式、正则型式的精确结果，找出符号矩阵的特征值而无由数值计算引入的误差。符号数学工具箱中的工具是建立在功能强大的称作Maple软件的基础上。它最初是由加拿大的滑铁卢（Waterloo）大学开发的。当要求MATLAB进行符号运算时，它就请求Maple去计算并将结果返回到MATLAB命令窗口。因此，在MATLAB中的符号运算是MATLAB处理数字的自然扩展。1.1符号计算基础一、定义符号变量参与符号运算的对象可以是符号变量、符号表达式或符号矩阵。符号变量要先定义，后引用。1.sym函数：每次定义一个符号变量。x=sym('x')%创建变量x。x=sym('x','real')%创建实变量x。x=sym('x','unreal')去掉x的实数属性。a=sym('alpha')%创建变量a，打印为alpha。2.syms函数sysmsarg1arg2…argN。一次可定义多个符号变量，使用起来更方便。注意各变量间用空格隔开。二、符号表达式符号表达式的书写格式与数值表达式一样，由符号变量、函数、运算符组成，如：symsabcxf=a*x^2+b*x+c也可以用f=sym('a*x^2+b*x+c')把符号表达式ax2+bx+c赋给f，但此时没有定义变量abcx。三、生成符号函数有两种方法创建符号函数，一是用符号表达式，如上面f=a*x^2+b*x+c定义了符号函数f，从而可以用diff等符号函数对f进行运算。如diff(f)将返回结果2*a*x+b。一是用M文件，如定义抽样函数：functionz=sinc(x)ifisequal(x,sym(0))%数值不能直接参与符号运算，先转化为符号元素。z=1;elsez=sin(x)/x;end四、默认符号变量\n54当字符表达式中含有多于一个的变量，但其中只有一个变量是独立变量时，如果不告诉MATLAB哪一个变量是独立变量，MATLAB将基于以下规则选择：选择在字母顺序中最接近x的字母，如果有两个距x相等的字母，就选择在字母表中较后的那一个。可以利用find(f,1)找出符号函数f中的默认符号变量。如：symsabxf=sin(a*x+b);findsym(f,1)%返回ans=x一、符号矩阵符号矩阵的元素为符号表达式，可用函数sym来产生符号矩阵。例：>>G=sym('[cos(t)，sin(t)；-sin(t)，cos(t)]')G=[cos(t)，sin(t)][-sin(t)，cos(t)]由于数值不能直接参与符号运算，因此先要将数值矩阵转化为符号矩阵。函数sym也可以把数值矩阵转换成符号形式。如：G=sym('[2/3sin(2);0.331/0.7]')G=[2/3,sin(2)][0.33,1/0.7]二、简化表达式1.simplify：Simplify是功能强大、通用的工具。它利用各种类型代数恒等式，包括求和、积分和分数幂、三角、指数和log函数、Bessel函数、超几何函数和函数，来简化表达式。例：symsxS=(x^2+5*x+6)/(x+2);R=simplify(S)R=x+32.expand：将表达式按幂由低到高展开成多项式形式（也可以是三角、指数、对数形式）。例：expand((x-2)*(x-4))ans=x^2-6*x+83.factor：因式分解，如factor(x^3-y^3)返回(x-y)*(x^2+x*y+y^2)4.collect：合并同类项，例：collect((exp(x)+x)*(x+2))ans=x^2+(exp(x)+2)*x+2*exp(x)5.numden：如果表达式是一个有理分式(两个多项式之比)，或是可以展开为有理分式(包括哪些分母为1的分式)，可利用numden来提取分子或分母6.horner(s)：嵌套，例：horner(x^3-6*x^2+11*x-6)ans=-6+(11+(-6+x)*x)*x7.simple：simple试用上面几种不同的简化方法，然后选择在结果表达式中含有最少字符的那种形式。三、可变精度运算符号工具箱支持可变精度运算，即支持符号计算并能以指定的精度返回结果。digits(5)%设置精度vpa('pi')%返回3.1416\n54一、函数计算器函数funtool格式funtool%该命令将生成三个图形窗口，FigureNo.1用于显示函数f的图形，FigureNo.2用于显示函数g的图形，FigureNo.3为一可视化的、可操作与显示一元函数的计算器界面。在该界面上由许多按钮，可以显示两个由用户输入的函数的计算结果：加、乘、微分等。funtool还有一函数存储器，允许用户将函数存入，以便后面调用。在开始时，funtool显示两个函数f(x)=x与g(x)=1在区间[-2*pi,2*pi]上的图形。Funtool同时在下面显示一控制面板，允许用户对函数f、g进行保存、更正、重新输入、联合与转换等操作。图8 1函数计算器1.2微积分一、极限函数limit格式limit(F,x,a)%计算符号表达式F=F(x)当x→a时的极限值。limit(F,a)%用命令findsym(F)确定F中的自变量，设为变量x，再计算F的极限值，当x→a时。limit(F)%用命令findsym(F)确定F中的自变量，设为变量x，再计算F的极限值，当x→0时。limit(F,x,a,'right')或limit(F,x,a,'left')%计算符号函数F的单侧极限：左极限x→a-或右极限x→a+。例：symsxath;Thenlimit(sin(x)/x)=>1limit(1/x,x,0,'right')=>inflimit((sin(x+h)-sin(x))/h,h,0)=>cos(x)二、导数（包括偏导数）函数diff格式diff(S,'v')、diff(S,sym('v'))%对表达式S中指定符号变量v计算S的1阶导数。diff(S)%对表达式S中的符号变量v计算S的1阶导数，其中v=findsym(S)。diff(S,n)%对表达式S中的符号变量v计算S的n阶导数，其中v=findsym(S)。diff(S,'v',n)%对表达式S中指定的符号变量v计算S的n阶导数。例：symsxtdiff(sin(x^2))%结果为2*cos(x^2)*x\n54diff(t^6,'t',6)%结果为returns720一、符号函数的积分函数int格式R=int(S,v)%对符号表达式S中指定的符号变量v计算不定积分。注意的是，表达式R只是函数S的一个原函数，后面没有带任意常数C。R=int(S)%对符号表达式S中的符号变量v计算不定积分，其中v=findsym(S)。R=int(S,v,a,b)%对表达式s中指定的符号变量v计算从a到b的定积分R=int(S,a,b)%对符号表达式s中的符号变量v计算从a到b的定积分，其中v=findsym(S)。例：symsxnint(x^n)%或int(x^n,x)、int(x^n,'x')，结果均为x^(n+1)/(n+1)int(sin(2*x),0,pi/2)%结果为1二、常微分方程的符号解函数dsolve格式r=dsolve('eq1,eq2,…','cond1,cond2,…','v')说明对给定的常微分方程（组）eq1,eq2,…中指定的符号自变量v，与给定的边界条件和初始条件cond1,cond2,….求符号解（即解析解）r；若没有指定变量v，则缺省变量为t；在微分方程（组）的表达式eq中，大写字母D表示对自变量(设为x)的微分算子：D=d/dx，D2=d2/dx2，…。微分算子D后面的字母则表示为因变量，即待求解的未知函数。初始和边界条件由字符串表示：y(a)=b，Dy(c)=d，D2y(e)=f，等等，分别表示，，；若边界条件少于方程（组）的阶数，则返回的结果r中会出现任意常数C1，C2，…；dsolve命令最多可以接受12个输入参量（包括方程组与定解条件个数，当然我们可以做到输入的方程个数多于12个，只要将多个方程置于一字符串内即可）。若没有给定输出参量，则在命令窗口显示解列表。若该命令找不到解析解，则返回一警告信息，同时返回一空的sym对象。这时，用户可以用命令ode23或ode45求解方程组的数值解。函数dsovle计算常微分方程的符号解。因为要求解微分方程，就需要用一种方法将微分包含在表达式中。所以，dsovle句法与大多数其它函数有一些不同，用字母D来表示求微分，D2，D3等等表示重复求微分，并以此来设定方程。任何D后所跟的字母为因变量。方程=0用符号表达式D2y=0来表示。独立变量可以指定或由symvar规则选定为缺省。例如，一阶方程dy/dx=1+y2的通解为：>>dsolve('Dy=1+y^2')%findthegeneralsolutionans=-tan(-x+C1)二阶方程的解dsolve('D2x+2*Dx+x=exp(-t)','x(0)=0','Dx(0)=0','t')ans=1/2*t^2*exp(-t)解微分方程组>>[f，g]=dsolve('Df=3*f+4*g'，'Dg=-4*f+3*g')f=C1*exp(3*x)*sin(4*x)+C2*exp(3*x)*cos(4*x)g=-C2*exp(3*x)*sin(4*x)+C1*exp(3*x)*cos(4*x)加上初始条件：f(0)=0和g(0)=1，我们可以得到：>>[f，g]=dsolve('Df=3*f+4*g'，'Dg=-4*f+3*g'，'f(0)=0，g(0)=1')f=exp(3*x)*sin(4*x)\n54g=exp(3*x)*cos(4*x)一、级数求和r=symsum(s,v)，求表达式s的级数和，变量v=0~v-1，缺省时变量为k例：symsknxsymsum(k^2)%return1/3*k^3-1/2*k^2+1/6*ksymsum(k)%returns1/2*k^2-1/2*k二、Taylor级数解析函数f(x)在点x=a的Taylor级数定义为：命令符号函数的Taylor级数展开式函数taylor格式r=taylor(f,n,v)%返回符号表达式f中的、指定的符号自变量v（若表达式f中有多个变量时）的n-1阶的Maclaurin多项式（即在零点附近v=0）近似式，其中v可以是字符串或符号变量。r=taylor(f)%返回符号表达式f中的、符号变量v的6阶的Maclaurin多项式（即在零点附近v=0）近似式，其中v=findsym(f)。r=taylor(f,n,v,a)%返回符号表达式f中的、指定的符号自变量v的n-1阶的Taylor级数（在指定的a点附近v=a）的展开式。其中a可以是一数值、符号、代表一数字值的字符串或未知变量。我们指出的是，用户可以以任意的次序输入参量n、v与a，命令taylor能从它们的位置与类型确定它们的目的。例：taylor(log(x),6,1)ans=x-1-1/2*(x-1)^2+1/3*(x-1)^3-1/4*(x-1)^4+1/5*(x-1)^51.2线性代数一、线性代数运算Matlab5.x版开始，“四则”符号运算指令形式上与数值计算完全相同，它们分别是+,-,*,/。加法：+（旧函数为symadd，仍可用）减法：-（旧函数为symsub）乘法：*（旧函数为symmul）除法：/（旧函数为symdiv）如：symsxx+3*x%返回ans=4*x二、矩阵函数建立对角矩阵（向量）：diag(v,k)上三角矩阵：triu(x,,k)下三角矩阵：tril(x,k)矩阵的逆阵：inv(A)矩阵的行列式：det(A)例：symsabcd;det([a,b;c,d])ans=a*d-b*c矩阵的秩：rank(A)\n54矩阵指数：expm(A)例：求微分方程组X'=AX的基解矩阵>>symstA=[21;-14];expm(A*t)ans=[exp(3*t)-t*exp(3*t),t*exp(3*t)][-t*exp(3*t),exp(3*t)+t*exp(3*t)]一、特征根和特征向量用函数eig可以求得符号矩阵的特征根和特征向量，格式为：lambda=eig(A)[V,D]=eig(A)二、矩阵的约当标准型函数jordan有两种格式jordan(F)%返回F的约当标准型矩阵[V，J]=jordan(F)%返回F的约当标准型矩阵J和相似变换矩阵V，使V-1FV=J1.2方程求解函数：g=solve(eq1,eq2,...,eqn,var1,var2,...,varn)功能：求解单个代数方程或代数方程组。如果表达式不是一个方程式(不含等号)，则在求解之前函数solve将表达式置成等于0。如果想对非缺省x变量求解，solve必须指定变量。>>solve('a*x^2+b*x+c'，'x')%solvefortherootsofthequadraticequtionans=[1/2/a*(-b+(b^2-4*a*c)^1/2)][1/2/a*(-b-(b^2-4*a*c)^1/2)]带有等号的符号方程也可以求解。>>f=solve('cos(x)=sin(x)')%solveforxf=1/4*pi用numeric可得到数值解。>>numeric(f)%返回ans=0.7854注意在求解周期函数方程时，有无穷多的解。在这种情况下，solve对解的搜索范围限制在接近于零的有限范围，并返回非唯一的解的子集。如果不能求得符号解，就计算可变精度解。>>x=solve('exp(x)=tan(x)')x=1.306326940423079代数方程组求解的例子：[x,y]=solve('x^2+x*y+y=3','x^2-4*x+3=0')x=[1][3]y=[1]\n54[-3/2]1.1积分变换一、傅立叶变换F=fourier(f,v)傅立叶变换f=ifourier(F,u)傅立叶逆变换例：symsxf=exp(-x^2);fourier(f)ans=pi^(1/2)*exp(-1/4*w^2)二、拉普拉斯变换拉普拉斯变换定义为：，因此拉普拉斯变换可看作频域中的傅里叶变换在复频域中的推广。若公式中的积分下限是0，则称作单边拉普拉斯变换。符号工具箱所求的是单边拉普拉斯变换。laplace(F,t)拉普拉斯变换F=ilaplace(L,y)拉普拉斯逆变换例：解微分方程x''+2x′+x=e-t，x(0)=x'(0)=0symstsXlaplace(exp(-t))%求e-t的拉普拉斯变换，ans=1/(s+1)X=solve('s*s*X+2*s*X+X-1/(s+1)','X')xt=ilaplace(X)%求X的拉普拉斯逆变换结果为X=1/(s+1)/(s^2+2*s+1)xt=1/2*t^2*exp(-t)%方程的解为三、Z变换Z变换用于研究离散时间系统，其定义为：，它可由采样信号的拉普拉斯变换推导出来。F=ztrans(f,z)Z变换f=iztrans(F,k)逆Z变换
