1、南京航空航天大学 王正盛1MATLAB 数学工具软件实例简明教程王 正 盛 编写南 京 航 空 航 天 大 学南京航空航天大学 王正盛2第一章 MATLAB简介MALAB译于矩阵实验室MATrix LABoratory是用来提供通往LINPACK和EISPACK矩阵软件包接口的 后来它渐渐发展成了通用科技计算图视交互系统和程序语言MATLAB的基本数据单位是矩阵它的指令表达与数学工程中常用的习惯形式十分相似比如矩阵方程Ax=b在MATLAB中被写成A*x=b而若要通过A,b求x那么只要写x=Ab即可完全不需要对矩阵的乘法和求逆进行编程因此用MATLAB解算问题要比用C Fortran等语言简捷
2、得多MATLAB发展到现在已经成为一个系列产品 MATLAB主包和各种可选的toolbox工具包主包中有数百个核心内部函数迄今所有的三十几个工具包又可分为两类功能性工具包和学科性工具包功能性工具包主要用来扩充MATLAB的符号计算功能图视建模仿真功能文字处理功能以及硬件实时交互功能这种功能性工具包用于多种学科而学科性工具包是专业性比较强的如控制工具包Control Toolbox信号处理工具包(SignalProcessing Toolbox) 通信工具包(Communication Toolbox)等都属此类开放性也许是MATLAB最重要最受人欢迎的特点除内部函数外所有MATLAB主包文件和
3、各工具包文件都是可读可改的源文件用户可通过对源文件的修改或加入自己编写文件去构成新的专用工具包MATLAB已经受了用户的多年考验在欧美发达国家MATLAB 已经成为应用线性代数自动控制理论数理统计数字信号处理时间序列分析动态系统仿真等高级课程的基本教学工具成为攻读学位的大学生硕士生博士生必须掌握的基本技能在设计研究单位和工业部门MATLAB被广泛地用于研究和解决各种具体工程问题第二章 MATLAB入门2.1 工作窗和指令行的操作除了通过菜单选项对工作窗进行控制外 MATLAB还提供了许多通过键盘输入的控制指令如下表MATLAB工作窗中的部分通用指令quit关闭和退出MATLABclc擦除MAT
4、LAB 工作窗中的所有显示内容clf擦除MATLAB 的当前图形窗中的图形clear清除内存中的变量和函数pack收集内存碎片以扩大内存空间dir列出指定目录下的文件和子目录清单cd改变当前工作子目录disp在运行中显示变量和文字内容type显示所有指定文件的全部内容echo控制运行文件指令是否显示的开关南京航空航天大学 王正盛3hold控制当前图形窗对象是否被刷新启动MATLAB 后就可以利用它工作了由于MATLAB 是一种交互式语言随时输入指令即时给出运算结果是它的主要工作方式当然更可以编制程序详见第七章比如要计算51302+).sin( 的值只要在光标位置处键入2*sin(0.3*pi)
5、/(1+sqrt(5) 然后按Enter键,该指令便被执行并给出结果ans = 0.5000 下面介绍控制光标对指令进行编辑的一些常用操作键常用操作键键 名作 用键 名作 用前寻式调回已输入过的指令行Home使光标移到当前行的首端后寻式调回已输入过的指令行End使光标移到当前行的尾端在 当前行中左移光标Delete删除光标右表边的字符在 当前行中右移光标Backspace删除光标左表边的字符PageUp前寻式翻阅当前窗中的内容Esc清楚当前行的全部内容PageDown后寻式翻阅当前窗中的内容2.2 简单矩阵的输入在MATLAB中矩阵输入的方法有多种此处只简单介绍矩阵的直接输入法详细介绍见第三章
6、在MATLAB中不必对矩阵维数做任何说明存储将自动配置在直接输入矩阵时矩阵元素用空格或逗号分隔矩阵行用隔离整个矩阵放在方括号 里例1A=1,2,3;4,5,6;7,8,9;10,11,12 A =1 2 34 5 67 8 910 11 12 说明指令执行后矩阵A被保存在MATLAB的工作间Workspace中以备后用如果用户不用clear指令清除它或对它重新定义该矩阵会一直保存在工作间中直到本MATLAB指令窗被关闭为止例2矩阵分行输入A=1 2 3 45 6 7 80 1 2 3 A = 1 2 3 45 6 7 8南京航空航天大学 王正盛40 1 2 3 例3矩阵元素输入B(1,2)=3
7、;B(4,4)=6;B(4,2)=11 B = 0 3 0 00 0 0 00 0 0 00 11 0 6 2.3 语句与变量MATLAB采用表达式语句用户输入语句由MATLAB系统结实运行MATLAB语句有两种常见的形式1表达式2变量=表达式说明1表达式由算符函数变量名和数字构成2在第一种形式中表达式被执行后产生的矩阵将被自动赋给名为ans的变量并 显示在屏幕上ans是一个缺省变量名它会被以后类似的操作刷新3在第二种形式中等号右边的表达式是被演绎后产生的矩阵将被赋给等号左边的变量存入内存并显示在屏幕上4 书写表达式时运算符号= +以及*等两侧允许有空格以增加可读性但在复数或符号表达式中要尽量
8、避免装饰性空格以防出错5变量名函数名以一个字母打头后面最多可接19个字母或数字注意MATLAB是区分字母的大小写的例1 表达式的计算结果2001/81 ans = 24.7037 例2 运算结果的赋值s=1-1/2+1/3-1/4+1/5-1/6+1/7-1/8; 说明结尾的分号作用是指令执行结果将不会显示在屏幕上但变量s仍将驻留在内存中如想看s的值只要键入s s = 0.6345 2.4 Who Whos和永久变量Who和Whos这两个指令的作用都是列出在MATLAB工作间中已经驻留的变量名清单不过Whos在给出变量名的同时还给出它们的维数及性质例1 用 who检查内存变量who Your
9、variables are:s 例2 用whos检查驻留变量的详细情况whos Name Size Bytes Classs 1x1 8 double arrayGrand total is 1 elements using 8 bytes 南京航空航天大学 王正盛5在MATLAB工作内存中还驻留几个由系统本身在启动时定义的变量如下表称为永久变量Permanent variables或称为预定义变量Predefined variables系统预定义的变量eps计算机的最小正数在pc机上它等于522pi圆周率的近似值3.14159265358979inf或Inf无穷大NaN不定量i,j虚数单位定
10、义1= jiflops浮点运算次数用于统计计算量说明1它们是在MATLAB启动时自定义的2它们不会被清除内存变量指令clear所清除 3他们可以重新定义为其他值但用clear可清除重定义值恢复预定义值例1 无穷大s=1/0 Warning: Divide by zero.s = Inf 无穷大a=Inf/inf a =NaN 2.5 数与表达式MATLAB的数值采用习惯的十进制表示可以带小数点或负号如下是合法的3 -99 0.0013 9.2445154 1.2434e-6 4.673e33在采用IEEE浮点算法的计算机上数值的相对精度是eps即大约保持16位有效数字数值范围大致为308308
11、101101 表达式由下列算符构成并按习惯的优先次序进行运算+ 加法 减法 * 乘法 / 右除 左除 乘方注意设置两种除法是为了方便矩阵的运算对标量而言两者作用相同例1x=2*pi/3+23/5-0.3e-3 x =3.6941 2.6 复数和复矩阵MATLAB认识复数并定义i和j作为虚数单位矩阵元素允许是复数复变量和由它们组成的表达式南京航空航天大学 王正盛6例1z1=3+4*i,z2=2*exp(i*pi/6)z=z1*z2 z1 =3.0000 + 4.0000iz2 =1.7321 + 1.0000iz =1.1962 + 9.9282i 例2A=1,3;2,4-i*5,8;6,9B=
12、1+5*i,2+6*i;3+8*i,4+9*iC=A*B A =1.0000 - 5.0000i 3.0000 - 8.0000i2.0000 - 6.0000i 4.0000 - 9.0000iB =1.0000 + 5.0000i 2.0000 + 6.0000i3.0000 + 8.0000i 4.0000 + 9.0000iC =1.0e+002 *0.9900 1.1600 - 0.0900i1.1600 + 0.0900i 1.3700 2.7 函数MATLAB的强大功能可函数中略见一斑本质上讲分为三类1 内部函数2 系统附带各种工具包中的M 文件所提供的大量函数3 用户自己增加的
13、函数这一特点是其他许多软件平台无法比拟的MATLAB提供的通用数理类函数包括z 基本数学函数z 特殊函数z 基本矩阵函数z 特殊矩阵函数z 矩阵分解和分析函数z 数据分析函数z 微分方程求解z 多项式函数z 非线性方程及其优化函数z 数值积分函数z 信号处理函数例z=1233.344x=sqrt(log(z) z =南京航空航天大学 王正盛71.233344000000000e+003x =2.66786140168028 2.8 显示格式在缺省的状态下MATLAB以短格式short格式显示计算结果可以用MATLAB命令窗口中format指令来改变数字的显示格式由于MATLAB以双精度执行所有
14、运算显示格式的设置仅影响矩阵的显示不影响矩阵的计算与存储如果矩阵的所有元素都是整数则矩阵以不带小数点的格式显示如果有一个元素不是整数则有几种输出格式默认格式为short格式只显示5位有效数字其他的显示格式可显示更多的有效数字还可用科学表示法例x=4/3 1.2345e-6默认short格式x =1.3333 0.0000 format short e 短格式科学表示x x =1.3333e+000 1.2345e-006 format long 长格式x x =1.33333333333333 0.00000123450000 format long e 长格式科学表示x x =1.33333
15、3333333333e+000 1.234500000000000e-006 format bank 银行格式x x = 1.33 0.00 format hex 十六进制格式x x = 3ff5555555555555 3eb4b6231abfd271 format + +格式用于显示大矩阵的紧凑格式+空格分别表示正数负数和零x x =+ 另外还有一种命令为format compact(紧凑格式)它消去了矩阵之间的间隔行这样可在一屏中显示更多的信息2.9 变量的存储与调用quit和exit指令都可退出MATLAB 结束MATLAB任务会删除工作间中的变南京航空航天大学 王正盛8量在退出前可以
16、保存工作空间以备再次调出使用这些变量保存的指令格式1 save 工作间中的所有变量保存在磁盘上名为matlab.mat的文件中2 save 文件名 变量名 将指定的变量保存在指定文件中如save temp x y z 把x,y,z这三个变量保存在文件temp.mat中在下次加载MATLAB时可以利用load指令将保存在文件中的变量恢复到工作间中其格式有1 load 将保存在matlab.mat中的变量装入到MATLAB 工作间中2 load 文件名 变量名 从指定的文件中将指定的变量装入MATLAB工作间如load temp x 从文件temp.mat中只将变量x装入到MATLAB工作间中2.
17、10 图形图形是MATLAB的主要特色之一MATLAB 图形指令具有自然简洁灵活及易扩充的特点MATLAB的指令很多这里仅介绍几个简单的绘图指令详见第六章例1 作多条曲线t=0:pi/50:4*pi;y0=exp(-t/3);y=exp(-t/3).*sin(3*t);plot(t,y,t,y0,t,-y0)grid0 2 4 6 8 10 12 14-1-0.8-0.6-0.4-0.200.20.40.60.81例2 三维曲面x=-8:0.5:8;y=x;X=ones(size(y)*x;Y=y*ones(size(x);R=sqrt(X.2+Y.2)+eps;Z=sin(R)./R;mes
18、h(Z);colormap(1,0,0) 南京航空航天大学 王正盛92.11 lp指令lookfor指令及其他帮助指令MATLAB的在线帮助系统相当完备就查询系统的调用方式而言可分为两种1 从MATLAB指令窗的 help菜单选项中寻求帮助此与一般windows的求助方法一样2 在MATLAB指令窗中直接键入求助指令(i)help 不带任何参数显示出MATLAB的目录项产生清单信息help HELP topics:matlabgeneral - General purpose commands.matlabops - Operators and special characters.matla
19、blang - Programming language constructs.matlabelmat - Elementary matrices and matrixmanipulation.matlabelfun - Elementary math functions.matlabspecfun - Specialized math functions.matlabmatfun - Matrix functions - numerical linearalgebra.matlabdatafun - Data analysis and Fouriertransforms.matlabpoly
20、fun - Interpolation and polynomials.matlabfunfun - Function functions and ODE solvers.matlabsparfun - Sparse matrices.matlabgraph2d - Two dimensional graphs.matlabgraph3d - Three dimensional graphs.matlabspecgraph - Specialized graphs.matlabgraphics - Handle Graphics.matlabuitools - Graphical user i
21、nterface tools.matlabstrfun - Character strings.matlabiofun - File input/output.matlabtimefun - Time and dates.matlabdatatypes - Data types and structures.matlabwinfun - Windows Operating System InterfaceFiles (DDE/ActiveX)matlabdemos - Examples and demonstrations.toolboxruntime - MATLAB Runtime Ser
22、ver DevelopmentKitrtwwindows - Real Time Windows Target.南京航空航天大学 王正盛10daqdaq - Data Acquisition Toolboxdaqdaqdemos - Data Acquisition Toolbox - DataAcquisition Demos.toolboxdials - Dials (2)在windows环境下建立MATLAB只能在启动时由mathabrc.m设定的路径上搜索不能与原定路径以外的其他目录交换信息可用以下三种方法扩充南京航空航天大学 王正盛14(1在MATLAB指令窗口中键入 CD C:MY
23、DIR2利用path指令扩展搜索路径 path(path, c:mydir)3在MATLAB 环境下键入pathtool 或者在MATLAB指令窗口菜单上File中的Set path项设置第三章 MATLAB的数值计算功能3.1 数值矩阵的创建保存和数据格式3.1.1 创建矩阵的直接输入法前面已述此出不赘述例x=14 ;y=4.32;A=x,2*x-y,0;sin(pi/4),3*y+x,sqrt(y) A =14.0000 23.6800 00.7071 26.9600 2.0785 3.1.2 利用MATLAB函数和语句创建数值矩阵 例1 利用指令reshape创建数值矩阵av=1:12
24、%产生12个元素的行向量av以%开头的是注释行bm=reshape(av,3,4) %利用向量av创建43矩阵bmav =1 2 3 4 5 6 7 8 9 10 1112bm =1 4 7 102 5 8 113 6 9 12 例2 利用指令diag产生对角阵ar=rand(4,4) %产生44的0-1均匀分布随即矩阵ard=diag(ar) %用矩阵的主对角线元素形成向量dD=diag(d) %用向量d构成对角矩阵Dar =0.9501 0.8913 0.8214 0.92180.2311 0.7621 0.4447 0.73820.6068 0.4565 0.6154 0.17630.4
25、860 0.0185 0.7919 0.4057d =0.95010.76210.61540.4057D =0.9501 0 0 00 0.7621 0 0南京航空航天大学 王正盛150 0 0.6154 00 0 0 0.4057 1.0.0 利用M文件创建和保存矩阵本节方法既适用于数值矩阵又适用于符号矩阵例1 创建和保存矩阵AM的matrix.m文件生成过程步骤1使用DOS的编辑器edit ,Windows的书写器(write)记事本notepad或其他字处理软件如Word等编辑如下AM=1 2 3;3 4 5步骤2把此内容以纯文本方式ASCII保存在用户自己的目录下名为matrix.m的
26、文件中步骤3在MATLAB指令窗中只要键入matrix矩阵AM 就会自动生成于MATLAB工作内存中即产生一个名为AM的变量供显示和调用3.1.4 通过MAT文件保存和获取矩阵MAT文件是MATLAB保存数据的一种标准格式二进制文件MAT文件的生成和调用由指令save和load进行例1 把矩阵AR保存到文件大他data.mat步骤1在矩阵AR 存在于MATLAB 内存空间的前提下键入save data AR步骤2在下次进入MATLAB后需要矩阵AR时键入如下边可将data.mat中的内容读入MATLAB 内存空间load data说明MATLAB默认扩展名为.mat默认路径为matlabbin
27、子目录用户如把data.mat登陆在指定目录可用如下命令保存或调入save c:mydirdata ARload c:mydirdata AR3.1.5 利用外部数据文件装入到指定矩阵假如磁盘中已有名为c:mydirdata.dat的ASCII数据文件利用指令load c:mydirdata.dat可在MATLAB工作空间产生一个名为 data的矩阵即变量当然也可以用指令fopen“ “fread“及其他MATLAB底层数据输入输出I/O指令实现可查看帮助如help fopen 3.2 矩阵的标识矩阵的元素子矩阵可以通过标识向量冒号的标识来援引和赋值例1b=1 2 3 4 5; 6 7 8 9
28、 10 ;11 12 13 14 15b23=b(2,3)b1=b(1:2,1 3 5)b2=b(3 1,:)b(1 3,2 4)=zeros(2) b =1 2 3 4 56 7 8 9 1011 12 13 14 15b23 =8南京航空航天大学 王正盛16b1 =1 3 56 8 10b2 =11 12 13 14 151 2 3 4 5b =1 0 3 0 56 7 8 9 1011 0 13 0 15 例2x=1 2 3 4 5 %产生51向量x =1 2 3 4 5l=xm时此方程称为超定方程Overdetermined Equation3) 当nm时此方程称为欠定方程Underd
29、etermined Equation3.5.1 矩阵逆和除法解恰定方程组1 采用求逆运算x=inv(A)*b2 采用左除运算x=Ab说明1由于MATLAB 遵循IEEE算法所以即使A阵奇异该运算也照样进行但在运算结束时一方面给出警告Warning:Matrix is singular to workingprecision另一方面所得逆阵的元素都是Inf无穷大1 当A为病态时也给出警告信息2 在MATLAB中inv指令不很有用MATLAB推荐尽量使用除运算少用逆运算例1求逆法和左除法解恰定方程组的性能对比为对比两种方法的性能先用以下指令构造一个条件数很大的高阶恰定方程组rand(seed,12
30、);%选定随机种子目的是可重复产生随机矩阵AA=rand(500)+1.e8;%rand(500)生成500500均匀分布的随机矩阵%每个随机矩阵元素加一个数的目的是使A的条件数变大x=ones(500,1); %令解向量x为全1的500元列向量b=A*x; %为使Ax=b方程一致用A和x生成b向量cond(A) %计算矩阵A的条件数ans =1.7608e+013过程1求逆法解恰定方程组的误差残差和所用计算时间tic %启动记时器Stopwatch Timerxi=inv(A)*b; %xi是用求逆法解恰定方程组所得的解toc %关闭计时器并显示解方程所用的时间eri=norm(x-xi)
31、%解向量xi与真解向量的2-范误差rei=norm(A*xi-b)/norm(b) %方程的2-范相对残差elapsed_time =7.2500eri =0.0066rei =1.5488e-006南京航空航天大学 王正盛20过程2左除法解恰定方程组的误差残差和所用计算时间ticxd=Ab;% xd是用左除法解恰定方程组所得的解tocerd=norm(x-xd)red=norm(A*xd-b)/norm(b) elapsed_time =3.3500erd =0.0021red =1.2695e-015 说明1 计算结果表明除法求解比求逆求解速度明显快精度相当但除法的相对残差几乎是机器零而逆
32、阵法的相对残差高得多2 MATLAB在设计求逆函数inv时采用的是Gauss消去法3 MATLAB在设计左除运算解恰定方程时并不求逆而是直接采用Gauss消去法求解有效地减少了残差所以即便在高条件数下也能得到较好的结果3.5.2 矩阵除法解超定方程组1 求正则方程Normal equations bAx)AA(TT=的解2 用Householder变换Householder transformation直接求原超定方程的最小二乘解由于第二种方程法采用的是正交变换据最小二乘理论可知第二种方法所得的解的准确性可靠性都比第一种方法好得多MATLAB解超定方程组用的就是第二种方法例除运算解超定方程的简
33、单算例a=1 2 3;4 5 -6;7 8 9;10 11 12;b=1:4;x=ab x = -0.33330.66670.00003.5.3 矩阵除法解欠定方程组欠定方程的解是不唯一的用除法运算所得的解有两个重要特征1在解中至多有Rank(A)个非零元素2它是这类型解中范数最小的一个例除运算解欠定方程的简单算例a=1 2 3;4 5 -6;7 8 9;10 11 12;b=a;c=1 3 3;x=bcx =2.00000.16670-0.1667南京航空航天大学 王正盛213.6 矩阵分解函数3.6.1 LU分解L,U=lu(X)产生一个上三角矩阵U和一个心理上下三角矩阵L即由下三角 矩阵
34、和置换矩阵构成使得X=L*U L,U,P=lu(X)产生一个上三角矩阵U和一个下三角矩阵L以及一个置换矩阵P使得P*X=L*U注意X必须是方阵例1用lu分解的两种指令格式对矩阵A进行 LU 分解解 A=1 1 1; 5 4 3; 2 1 1L,U=lu(A)L,U,P=lu(A) A =1 -1 15 -4 32 1 1L,U=lu(a)L =0.2000 -0.0769 1.00001.0000 0 00.4000 1.0000 0U =5.0000 -4.0000 3.00000 2.6000 -0.20000 0 0.3846L,U,P=lu(a)L =1.0000 0 00.4000
35、1.0000 00.2000 -0.0769 1.0000U =5.0000 -4.0000 3.00000 2.6000 -0.20000 0 0.3846P =0 1 00 0 11 0 0例2 A为严格对角占优矩阵A=4 1 2 ;2 5 1; 5 3 11L,U=lu(a) A =4 1 22 5 -15 3 11L,U=lu(A)南京航空航天大学 王正盛22L=0.8000 -0.3684 1.00000.4000 1.0000 01.0000 0 0U =5.0000 3.0000 11.00000 3.8000 -5.40000 0 -8.7895例3利用 LU 分解求解线性方程
36、组AX=bA=1,-1,1;5,-4,3;2,1,1 A = 1 -1 15 -4 32 1 1b=2;-3;1b = 2-31L,U,P=lu(A)L =1.0000 0 00.4000 1.0000 00.2000 -0.0769 1.0000U =5.0000 -4.0000 3.00000 2.6000 -0.20000 0 0.3846P =0 1 00 0 11 0 0y=L(P*b)y =-3.00002.20002.7692x=Uyx =-3.80001.40007.20001.0.0 QR分解Q,R=qr(A)产生一个与矩阵A相同维数的上三角矩阵R和一个酉阵 Q使得A=Q*R
37、Q,R,E=qr(A)产生一个置换矩阵E 上三角矩阵R和一个酉阵 Q使得A*E=Q*R例1a=1 1 1;5 4 3;2 1 1南京航空航天大学 王正盛23Q,R=qr(a)Q,R,E=qr(a) a =1 -1 15 -4 32 1 1Q,R=qr(a)Q = -0.1826 -0.1501 -0.9717-0.9129 -0.3412 0.2242-0.3651 0.9279 -0.0747R =-5.4772 3.4689 -3.28630 2.4427 -0.24560 0 -0.3737Q,R,E=qr(a)Q =-0.1826 -0.1501 -0.9717-0.9129 -0.3
38、412 0.2242-0.3651 0.9279 -0.0747R =-5.4772 3.4689 -3.28630 2.4427 -0.24560 0 -0.3737E =1 0 00 1 00 0 1例2 a为严格对角占优矩阵a=10 1 1;5 14 3;2 1 21Q,R=qr(a)Q =-0.1826 -0.1501 -0.9717-0.9129 -0.3412 0.2242-0.3651 0.9279 -0.0747R =-5.4772 3.4689 -3.28630 2.4427 -0.24560 0 -0.3737 Q,R,E=qr(a)Q =-0.0471 -0.0678 -
39、0.9966-0.1413 -0.9872 0.0738-0.9889 0.1443 0.0369R =-21.2368 1.0359 -3.15490 14.0331 -5.32540 0 -9.5230E =0 0 1南京航空航天大学 王正盛240 1 01 0 03.6.3 Choleshy分解Cholesky分解要求被分解矩阵A为对称正定矩阵R=chol(A)产生一个上三角矩阵R使得R*RAT=例a=2,1,1;1,2,-1;1,-1,3 R=chol(a)a =2 1 11 2 -11 -1 3R=chol(a)R =1.4142 0.7071 0.70710 1.2247 -1.2
40、2470 0 1.00003.7 多项式3.7.1多项式的表达和创建多项式表达方式的约定多项式nnnnaxaxaxa)x(p +=1110L 用以下系数行向量表示a,a,a,apnn 110 = L多项式行向量的创建方法1 多项式系数向量的直接输入法2 利用指令p=poly(AR)产生多项式系数向量若AR是方阵则多项式为特征多项式若AR是向量即ar,ar,arARnL21=则所得的多项式满足= )arx()arx)(arx(nL21 nnnnaxaxaxa +1110L例求3阶方阵A的特征多项式a=11 12 13;14 15 16;17 18 19;pa=poly(a)ppa=poly2st
41、r(pa,s) pa =1.0000 -45.0000 -18.0000 0.0000ppa =s3 - 45 s2 - 18 s + 5.3518e-015 注意1 n阶方阵的特征多项式系数向量一定是(n+1)维的2特征多项式向量的第一的元素必是1例由给定根向量求多项式系数向量r=-0.5,-0.3+0.4*i,-0.3-0.4*i;p=poly(r)pr=real(p)南京航空航天大学 王正盛25ppr=poly2str(pr,x) p = 1.0000 1.1000 0.5500 0.1250pr = 1.0000 1.1000 0.5500 0.1250ppr = x3 + 1.1 x
42、2 + 0.55 x + 0.125 说明1 要形成实数多项式根向量中的复数根必须共轭成对2 含复数根的根向量所生成的多项式系数向量如p的系数有可能带在截断误差数量级的虚部此时可采用取实部的指令real把虚部去掉3 poly2str是一个函数文件它存在于MATLAB的控制工具包(ControlToolbox)中3.6.2 常用多项式运算指令R=roots(p) 求多项式向量p的根PA=polyval(p,S) 按数组运算规则计算多项式值p为多项式S为矩阵PM=polyvalm(p;S) 按矩阵运算规则计算多项式值p为多项式S为矩阵r,p,k=residue(b,a) 部分分式展开b,a分别是分
43、子分母多项式系数向量r,p,k分别是留数极点直项向量P=polyfit(x,y,n) 用n阶多项式拟合x,y向量给定的数据例1求多项式2772623 xxx的根解r=roots(1,-6,-72,-27) r =12.1229-5.7345-0.3884 注意MATLAB约定多项式系数用行向量表示一组根用列向量表示例2两种多项式求值指令的差别s=pascal(4)p=poly(s);pp=poly2str(p,x)pa=polyval(p,s)pm=polyvalm(p,s) s =1 1 1 11 2 3 41 3 6 101 4 10 20pp =x4 - 29 x3 + 72 x2 -
44、29 x + 1南京航空航天大学 王正盛26pa =1.0e+004 *0.0016 0.0016 0.0016 0.00160.0016 0.0015 -0.0140 -0.05630.0016 -0.0140 -0.2549 -1.20890.0016 -0.0563 -1.2089 -4.3779pm =1.0e-011 *-0.0077 0.0053 -0.0096 0.0430-0.0068 0.0481 -0.0110 0.12220.0075 0.1400 -0.0095 0.26080.0430 0.2920 -0.0007 0.4737 说明: pm 中的元素都很小它是运算误
45、差造成的理论上它应该是零这就是著名的Caylay-Hamilton定理任何一个矩阵满足它自己的特征多项式例3用6阶多项式对区间0 2.5上的误差函数dte)x(yxt=022进行最小二乘拟合解x=0:0.1:2.5;y=erf(x); %计算误差函数在0,2.5内的数据点p=polyfit(x,y,6)px=poly2str(p,x)p =0.0084 -0.0983 0.4217 -0.7435 0.1471 1.10640.0004px =0.0084194 x6 - 0.0983 x5 + 0.42174 x4 - 0.74346 x3+ 0.1471 x2+ 1.1064 x + 0.
46、00044117 例4有效拟合的区间性图示用0,2.5区间数据拟合曲线拟合0,5区间数据x=0:0.1:5;x1=0:0.1:2.5y=erf(x);y1=erf(x1);p=polyfit(x1,y1,6)f=polyval(p,x);plot(x,y,bo,x,f,r-)axis(0,5,0,2)legend(拟合曲线,原数据线) x1 =Columns 1 through 70 0.1000 0.2000 0.3000 0.4000 0.50000.6000Columns 8 through 140.7000 0.8000 0.9000 1.0000 1.1000 1.20001.3000Columns 15 through 21南京航空航天大学
