matlab数值分析实验

数值分析实验董海云数理学院数学实验教学中心目录0Matlab介绍入门知识 31绪论 171.1例题解答 171.2Matlab中数值计算精度 202线性方程组的直接解法 222.1例题解答 222.2Matlab解线性方程组常用命令介绍 363线性方程组的迭代解法 383.1例题解答 383.2Matlab迭代解法用到的函数介绍 534方阵特征值和特征向量的计算 554.1例题解答 554.2Matlab关于方阵特征值为特征向量函数介绍 625非线性方程求根 645.1例题解答 645.2Matlab非线性方程求根的命令 856插值法 866.1例题解答 866.2Matlab插值函数介绍 1017数据拟合和最佳平方逼近 1037.1例题解答 1037.2Matlab数据拟合命令介绍 1138数值积分与数值微分 1148.1例题解答 1149常微分方程数值解法 1389.1例题解答 1389.2Matlab常微分方程数值解常用命令介绍 1540Matlab介绍入门知识1.Matlab简介MATLAB的含义是矩阵实验室（MATRIXLABORATORY），主要用于方便矩阵的存取，其基本元素是无须定义维数的矩阵.MATLAB自问世以来,就是以数值计算称.MATLAB进行数值计算的基本单位是复数数组（或称阵列），这使得MATLAB高度“向量化”.经过十几年的完善和扩充，现已发展成为线性代数课程的标准工具.由于它不需定义数组的维数，并给出矩阵函数、特殊矩阵专门的库函数，使之在求解诸如信号处理、建模、系统识别、控制、优化等领域的问题时，显得大为简捷、高效、方便，这是其它高级语言所不能比拟的.MATLAB中包括了被称作工具箱（TOOLBOX）的各类应用问题的求解工具.工具箱实际上是对MATLAB进行扩展应用的一系列MATLAB函数（称为M文件），它可用来求解各类学科的问题，包括信号处理、图象处理、控制系统辨识、神经网络等.随着MATLAB版本的不断升级，其所含的工具箱的功能也越来越丰富，因此，应用范围也越来越广泛，MATLAB提供的工具箱已覆盖信号处理、系统控制、统计计算、优化计算、神经网络、小波分析、偏微分方程、模糊逻辑、动态系统模拟、系统辨识和符号运算等领域.当前它的使用范围涵盖了工业、电子、医疗、建筑等各行各业.MATLAB中包括了图形界面编辑GUI，让使用者也可以象VB、VC、VJ、DELPHI等那样进行一般的可视化的程序编辑.在命令窗口（matlabcommandwindow）键入simulink，就出现(SIMULINK)窗口.以往十分困难的系统仿真问题，用SIMULINK只需拖动鼠标即可轻而易举地解决问题，这也是近来受到重视的原因所在.MATLAB语言由美国TheMathWorks开发,最早是由C.Moler用Fortran语言编写的，用来方便地调用LINPACK和EISPACK矩阵代数软件包的程序.后来他创立了MATHHWORKS公司，对MATLAB作了大量的、坚持不懈的改进.CleveB.Moler是TheMathWork公司的主席和首席科学家.曾任密歇系教授.他在两个计算机硬件制造商Intel公司的Hypercube组织和ArdenComputers公司工作了五年.他的主要专业兴趣在于数值分析和科学计算.他是MATLAB软件的创始者，也是著名的矩阵计算软件包LINPACK和EISPACK的著作这一，已撰写了三本有相关数值方法的教材.同时，他在SIAM（美国工业与应用数学学会）历任期刊编辑、委员会成员和副总裁，并从1996年开始担任理事会成员.2.Matlab入门知识Matlab变量名是以字母开头，后接字母、数字或下划线的字符序列，最多63个字符.在MATLAB中，变量名区分字母的大小写.赋值语句:变量=表达式或表达式其中表达式是用运算符将有关运算量连接起来的式子，其结果是一个矩阵.clear命令用于删除MATLAB工作空间中的变量.who和whos这两个命令用于显示在MATLAB工作空间中已经驻留的变量名清单.who命令只显示出驻留变量的名称，whos在给出变量名的同时，还给出它们的大小、所占字节数及数据类型等信息.利用MAT文件可以把当前MATLAB工作空间中的一些有用变量长久地保留下来，扩展名是.mat.MAT文件的生成和装入由save和load命令来完成.常用格式为：save文件名[变量名表][-append][-ascii]load文件名[变量名表][-ascii]其中，文件名可以带路径，但不需带扩展名.mat，命令隐含一定对.mat文件进行操作.变量名表中的变量个数不限，只要内存或文件中存在即可，变量名之间以空格分隔.当变量名表省略时，保存或装入全部变量.-ascii选项使文件以ASCII格式处理，省略该选项时文件将以二进制格式处理.save命令中的-append选项控制将变量追加到MAT文件中.(1)向量的创建用步长生成法：数组=初值:步长(增量):终值>>a=1:0.5:3a=1.00001.50002.00002.50003.0000用linspace生成：数组=linspace(初值,终值,等分点数目)>>b=linspace(1,3,5)b=1.00001.50002.00002.50003.0000列向量用分号（;）作为分行标记：>>c=[1;2;3;4;]c=1234若不想输出结果，在每一条语句后用分号作为结束符，若留空或用逗号结束，则在执行该语句后会把结果输出来.>>a+b;>>a+bans=23456(2)矩阵的创建直接输入:最简单的建立矩阵的方法是从键盘直接输入矩阵的元素.具体方法如下：将矩阵的元素用方括号括起来，按矩阵行的顺序输入各元素，同一行的各元素之间用空格或逗号分隔，不同行的元素之间用分号分隔.>>A=[123;456;235]A=123456235利用矩阵函数创建:>>B=magic(3)%魔方阵B=816357492>>C=hilb(3)%3阶Hilbert矩阵C=1.00000.50000.33330.50000.33330.25000.33330.25000.2000Matlab中用%引导注释其它创建矩阵函数还有:eye(m,n):生成m行n列单位矩阵.zeros(m,n):生成m行n列全零矩阵.ones(m,n):生成全1矩阵.rand(m,n):生成随机矩阵.rand:生成一个随机数.diag(A):取A的对角线元素.tril(A):取A的下三角元素.triu(A):取A的上三角元素.hilb(n):生成n维Hilbert矩阵.randn(n)：产生均值为0，方差为1的标准正态分布随机矩阵.vander(V):生成以向量V为基础向量的范得蒙矩阵.invhilb(n):求n阶的希尔伯特矩阵的逆矩阵.toeplitz(x,y):生成一个以x为第一列，y为第一行的托普利兹矩阵compan(p):生成伴随矩阵,p是一个多项式的系数向量，高次幂系数排在前，低次幂排在后.pascal(n):生成一个n阶帕斯卡矩阵.compan:生成伴随矩阵(3)矩阵运算MATLAB的基本算术运算有：＋(加)、－(减)、*(乘)、/(右除)、\(左除)、^(乘方).加法：>>A+Bans=939710136127减法：>>B-Aans=7-13-10126-3乘法：>>A*Bans=263826718371456243除法:>>magic(3)/hilb(3)ans=1.0e+003*0.2160-1.17601.14000.0570-0.40800.4500-0.22801.2240-1.1400在MATLAB中，有一种特殊的运算，因为其运算符是在有关算术运算符前面加点，所以叫点运算.点运算符有.*、./、.\和.^.两矩阵进行点运算是指它们的对应元素进行相关运算，要求两矩阵的维参数相同.>>A.*Bans=821812254282710MATLAB提供了6种关系运算符：<(小于)、<=(小于或等于)、>(大于)、>=(大于或等于)、==(等于)、～=(不等于).>>A>Bans=010100001MATLAB提供了3种逻辑运算符：&(与)、|(或)和～(非).在逻辑运算中，确认非零元素为真，用1表示，零元素为假，用0表示.设参与逻辑运算的是两个标量a和b，那么，a&ba,b全为非零时，运算结果为1，否则为0.a|ba,b中只要有一个非零，运算结果为1.～a当a是零时，运算结果为1；当a非零时，运算结果为0.3.矩阵操作和矩阵函数矩阵通过下标引用矩阵的元素，矩阵元素的序号就是相应元素在内存中的排列顺序.在MATLAB中，矩阵元素按列存储，先第一列，再第二列，依次类推.(1)矩阵拆分利用冒号表达式获得子矩阵.A(:,j)表示取A矩阵的第j列全部元素；A(i,:)表示A矩阵第i行的全部元素；A(i,j)表示取A矩阵第i行、第j列的元素.A(i:i+m,:)表示取A矩阵第i～i+m行的全部元素；A(:,k:k+m)表示取A矩阵第k～k+m列的全部元素，A(i:i+m,k:k+m)表示取A矩阵第i～i+m行内，并在第k～k+m列中的所有元素.此外，还可利用一般向量和end运算符来表示矩阵下标，从而获得子矩阵.end表示某一维的末尾元素下标.(2)利用空矩阵删除矩阵的元素在MATLAB中，定义[]为空矩阵.给变量X赋空矩阵的语句为X=[].(3)矩阵的转置转置运算符是单撇号(‘).(4)矩阵的旋转利用函数rot90(A,k)将矩阵A旋转90o的k倍，当k为1时可省略.(5)矩阵的左右翻转对矩阵实施左右翻转是将原矩阵的第一列和最后一列调换，第二列和倒数第二列调换，…，依次类推.MATLAB对矩阵A实施左右翻转的函数是fliplr(A).(6)矩阵的上下翻转MATLAB对矩阵A实施上下翻转的函数是flipud(A).(7)方阵A的逆矩阵inv(A)>>A=magic(3)A=816357492>>B=inv(A)B=0.1472-0.14440.0639-0.06110.02220.1056-0.01940.1889-0.1028>>A*Bans=1.00000-0.0000-0.00001.000000.000001.0000(8)方阵的行列式>>det(A)ans=-360(9)矩阵的迹>>C=trace(A)C=15(10)一些常用的基本初等三角函数三角函数:sin(x),cos(x),tan(x)反三角函数:asin(x),acos(x),atan(x)指数函数:exp(x)自然对数:log(x)常用对数:log10(x)以2为底的对数:log2(x)开平方:sqrt(x)绝对值:abs(x)计算一般函数值:eval(f)求虚部函数:imag(x)求实部函数:real(x)角相位函数:angle(x)共轭复数函数:conj(x)沿零方向取整:fix(x)舍入取整:round(x)沿负无穷大方向取整:floor(x)沿正无穷大方向取整:ceil(x)求除法的余数:rem符号函数:sign(x)最大公约数:gcd()4.图形可视化(1)二维绘图指令plotplot函数的基本调用格式为:plot(x,y,)其中x和y为长度相同的向量，分别用于存储x坐标和y坐标数据.plot(x)plot函数最简单的调用格式.当x是实向量时，以该向量元素的下标为横坐标，元素值为纵坐标画出一条连续曲线.实际上是绘制折线图.plot(x1,y1,x2,y2,…,xn,yn)当输入参数都为向量时，x1和y1，x2和y2，…，xn和yn分别组成一组向量对，每一组向量对的长度可以不同.每一向量对可以绘制出一条曲线，可以在同一坐标内绘制出多条曲线.plotyy(x1,y1,x2,y2)绘制出具有不同纵坐标标度的两个图形.holdon/off保持原有图形还是刷新原有图形,不带参数的hold命令在两种状态之间进行切换.plot(x1,y1,选项1,x2,y2,选项2,…,xn,yn,选项n)设置曲线样式进行绘图.选项字段见下表:表Matlab常用线形与颜色标记表符号线型符号线型符号线型颜色含义-实线.实点标记^朝上三角g绿色--虚线o圆圈<朝左三角r红色：点线X叉字符>朝右三角c青色-.点划线+加号p五角星m洋红*星号h六角形y黄色s方块k黑色d菱形w白色v朝下三角b蓝色(2)图形标注:title('图形名称'):图形标题xlabel('x轴说明')ylabel('y轴说明')text(x,y,'图形说明')legend('图例1','图例2',…)gtext('用鼠标确定位置的字符说明')(3)坐标控制axisaxis([xminxmaxyminymaxzminzmax])axis函数功能丰富，常用的格式还有：axisequal：纵、横坐标轴采用等长刻度.axissquare：产生正方形坐标系(缺省为矩形).axisauto：使用缺省设置.axisoff：取消坐标轴.axison：显示坐标轴.gridon/off:网格开/关boxon/off:加/不加边框线上述命令示例如下:>>x=1:length(peaks);>>plot(x,peaks);>>boxon;>>title('绘制混合图形');>>xlabel('X轴');>>ylabel('Y轴');绘制图像为:(4)二维数值函数的专用绘图函数fplotfplot(functionname,[a,b],tol,选项)其中functionname为函数名，以字符串形式出现，[a,b]为绘图区间，tol为相对允许误差，其系统默认值为2e-3.选项定义与plot函数相同.>>fplot(@(x)[tan(x),sin(x),cos(x)],2*pi*[-11-11]);(5)二维符号函数曲线专用命令ezplotf=f(x)时:ezplot(f)：在默认区间-2π<x<2π绘制f=f(x)的图形.ezplot(f,[a,b])：在区间a<x<b绘制f=f(x)的图形f=f(x,y)时:ezplot(f)：在默认区间-2π<x<2π和-2π<y<2π绘制f(x,y)=0的图形.ezplot(f,[xmin,xmax,ymin,ymax])：在区间xmin<x<xmax和ymin<y<ymax绘制f(x,y)=0的图形ezplot(f,[a,b])：在区间a<x<b和a<y<b绘制f(x,y)=0的图形若x=x(t),y=y(t):ezplot(x,y)：在默认区间0<t<2π绘制x=x(t)和y=y(t)的图形.ezplot(x,y,[tmin,tmax])：在区间tmin<t<tmax绘制x=x(t)和y=y(t)的图形>>figure;ezplot('cos(tan(pi*x))',[0,1]);(6)图形窗口的分割subplotsubplot(m,n,p)该函数将当前图形窗口分成m×n个绘图区，即每行n个，共m行，区号按行优先编号，且选定第p个区为当前活动区.在每一个绘图区允许以不同的坐标系单独绘制图形.(7)其他坐标系下的二维数据曲线图对数坐标图形:semilogx(x1,y1,选项1,x2,y2,选项2,…)semilogy(x1,y1,选项1,x2,y2,选项2,…)loglog(x1,y1,选项1,x2,y2,选项2,…)极坐标图polar:polar(theta,r,选项)其中theta为极坐标极角，r为极坐标矢径，选项的内容与plot函数相似.二维统计分析图:bar(x,y,选项):条形图stairs(x,y,选项):阶梯图stem(x,y,选项):杆图fill(x1,y1,选项1,x2,y2,选项2,…):填充图(8)三维曲线plot3plot3(x1,y1,z1,选项1,x2,y2,z2,选项2,…,xn,yn,zn,选项n)其中每一组x,y,z组成一组曲线的坐标参数，选项的定义和plot函数相同.当x,y,z是同维向量时，则x,y,z对应元素构成一条三维曲线.当x,y,z是同维矩阵时，则以x,y,z对应列元素绘制三维曲线，曲线条数等于矩阵列数.>>t=0:0.1:8*pi;>>plot3(sin(t),cos(t),t);(9)产生三维数据在MATLAB中，利用meshgrid函数产生平面区域内的网格坐标矩阵.其格式为：[X,Y]=meshgrid(x,y);语句执行后，矩阵X的每一行都是向量x，行数等于向量y的元素的个数，矩阵Y的每一列都是向量y，列数等于向量x的元素的个数.(10)绘制三维曲面的函数surf函数和mesh函数的调用格式为:mesh(x,y,z,c)surf(x,y,z,c)一般情况下，x,y,z是维数相同的矩阵.x,y是网格坐标矩阵，z是网格点上的高度矩阵，c用于指定在不同高度下的颜色范围.(11)标准三维曲面sphere函数的调用格式为：[x,y,z]=sphere(n)cylinder函数的调用格式为：[x,y,z]=cylinder(R,n)MATLAB还有一个peaks函数，称为多峰函数，常用于三维曲面的演示.(12)其他三维绘图指令介绍bar3函数绘制三维条形图，常用格式为bar3(y)bar3(x,y)stem3函数绘制离散序列数据的三维杆图，常用格式为：stem3(z)stem3(x,y,z)pie3函数绘制三维饼图，常用格式为：pie3(x)fill3函数等效于三维函数fill，可在三维空间内绘制出填充过的多边形，常用格式为：fill3(x,y,z,c)5.程序控制结构(1)数据的输入:A=input(提示信息，选项)其中提示信息为一个字符串，用于提示用户输入什么样的数据.如果在input函数调用时采用's'选项，则允许用户输入一个字符串.(2)数据的输出:disp(输出项)(3)程序的暂停:pause(延迟秒数)如果省略延迟时间，直接使用pause，则将暂停程序，直到用户按任一键后程序继续执行.若要强行中止程序的运行可使用Ctrl+C命令.(4)单分支if语句：if条件语句组end当条件成立时，则执行语句组，执行完之后继续执行if语句的后继语句，若条件不成立，则直接执行if语句的后继语句.(5)双分支if语句：if条件语句组1else语句组2end当条件成立时，执行语句组1，否则执行语句组2，语句组1或语句组2执行后，再执行if语句的后继语句.(6)多分支if语句：if条件1语句组1elseif条件2语句组2……elseif条件m语句组melse语句组nend语句用于实现多分支选择结构.(7)switch语句:switch表达式case表达式1语句组1case表达式2语句组2……case表达式m语句组motherwise语句组nend(8)try语句语句格式为：try语句组1catch语句组2endtry语句先试探性执行语句组1，如果语句组1在执行过程中出现错误，则将错误信息赋给保留的lasterr变量，并转去执行语句组2.(9)for语句for语句的格式为：for循环变量=表达式1:表达式2:表达式3循环体语句end其中表达式1的值为循环变量的初值，表达式2的值为步长，表达式3的值为循环变量的终值.步长为1时，表达式2可以省略.for语句更一般的格式为：for循环变量=矩阵表达式循环体语句end执行过程是依次将矩阵的各列元素赋给循环变量，然后执行循环体语句，直至各列元素处理完毕.(10)while语句while语句的一般格式为：while(条件)循环体语句end其执行过程为：若条件成立，则执行循环体语句，执行后再判断条件是否成立，如果不成立则跳出循环.(11)break语句和continue语句与循环结构相关的语句还有break语句和continue语句.它们一般与if语句配合使用.break语句用于终止循环的执行.当在循环体内执行到该语句时，程序将跳出循环，继续执行循环语句的下一语句.continue语句控制跳过循环体中的某些语句.当在循环体内执行到该语句时，程序将跳过循环体中所有剩下的语句，继续下一次循环.(12)循环的嵌套如果一个循环结构的循环体又包括一个循环结构，就称为循环的嵌套，或称为多重循环结构.(13)函数文件的基本结构函数文件由function语句引导，其基本结构为function输出形参表=函数名(输入形参表)注释说明部分函数体语句其中以function开头的一行为引导行，表示该M文件是一个函数文件.函数名的命名规则与变量名相同.输入形参为函数的输入参数，输出形参为函数的输出参数.当输出形参多于一个时，则应该用方括号括起来.(14)函数调用函数调用的一般格式是：[输出实参表]=函数名(输入实参表)注意的是，函数调用时各实参出现的顺序、个数，应与函数定义时形参的顺序、个数一致，否则会出错.函数调用时，先将实参传递给相应的形参，从而实现参数传递，然后再执行函数的功能.在MATLAB中，函数可以嵌套调用，即一个函数可以调用别的函数，甚至调用它自身.一个函数调用它自身称为函数的递归调用.(15)函数参数的可调性在调用函数时，MATLAB用两个永久变量nargin和nargout分别记录调用该函数时的输入实参和输出实参的个数.只要在函数文件中包含这两个变量，就可以准确地知道该函数文件被调用时的输入输出参数个数，从而决定函数如何进行处理.(16)全局变量与局部变量全局变量用global命令定义，格式为：global变量名(17)程序调试Debug菜单项:该菜单项用于程序调试，需要与Breakpoints菜单项配合使用.Breakpoints菜单项:该菜单项共有6个菜单命令，前两个是用于在程序中设置和清除断点的，后4个是设置停止条件的，用于临时停止M文件的执行，并给用户一个检查局部变量的机会，相当于在M文件指定的行号前加入了一个keyboard命令.调试命令:除了采用调试器调试程序外，MATLAB还提供了一些命令用于程序调试.命令的功能和调试器菜单命令类似，具体使用方法请读者查询MATLAB帮助文档.1绪论1.1例题解答例1.1计算,.解:创建符号函数:>>symsx;>>f=sym('sin(x)')f=sin(x)展开至7阶泰勒级数：>>h=taylor(f,8,0)h=x-1/6*x^3+1/120*x^5-1/5040*x^7求泰勒级数在处的函数值：>>subs(h,x,0.5)ans=0.479425533234127也可以通过内联函数来求解:>>H=inline(h)H=Inlinefunction:H(x)=x-1./6.*x.^3+1./120.*x.^5-1./5040.*x.^7>>feval(H,0.5)ans=0.479425533234127例1.2计算积分值.解:解法一:(符号法)：>>I=int('1/(1+x)','x',0,1)I=log(2)解法二:(数值法)：>>x=0:0.2:1;%将[0,1]等分为4等份>>f=1./(1+x);%分别计算每一个等分点的函数值>>I=0;>>fori=1:5I=I+(f(i)+f(i+1))/2*0.2;%将每一小曲边的梯形累加起来作为积分值End>>vpa(I,9)%取结果的小数精度为9位小数ans=.695634921例1.3略例1.4不用开平方根计算的值.解:解法一(符号法):>>A=sym('a');>>sqrt(A)ans=a^(1/2)解法二(数值法):按以下迭代公式迭代计算近似值:建立函数文件msqrt.mfunctionx=msqrt(x0,a)%用迭代法近似计算平方根%x0为初始迭代值,a为开平方数formatlong;x=zeros(20,1);x(1)=x0;fori=2:20x(i)=1/2*(x(i-1)+a/x(i-1));enddisp(x);用编写的函数计算,:>>msqrt(2,3);2.0000000000000001.7500000000000001.7320508100147271.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.732050807568877上述结果为迭代过程计算的中间结果,分析数据可知迭代收敛速度快,只需四次计算即可计算出较为准确的数值.例1.5略.例1.6计算,视已知数为精确数,用4位浮点数计算.解:直接在Matlab中输入式子:>>1/759-1/760ans=1.7336e-006若先转化为浮点数再运算可得:>>a=1/759,b=1/760,a-ba=0.0013b=0.0013ans=1.7336e-006可见Matlba在计算时,数据结构都取为双精度而提高了运算准确度.若以符号运算计算之,有:>>a=sym('1/759'),b=sym('1/760'),c=a-ba=1/759b=1/760c=1/576840可见符号运算准确但耗费运算时间.例1.7略.例1.8解方程.解:符号法解方程:>>x=solve('x^2-18*x+1','x')x=9+4*5^(1/2)9-4*5^(1/2)将结果保留小数点6位:>>vpa(x,6)ans=17.9443.5572e-11.2Matlab中数值计算精度1.Matlab中有三种运算精度,它们分别为数值算法、符号算法和可控精度算法,将它们分别介绍如下:数值算法把每个数取为16位,计算按浮点运算进行,它是运算速度最快的一种算法.符号算法把每个数都变为符号量,运算按有理量计算进行,它的优点是能够得到精确结果,缺点是占用空间大,并且运算速度最慢.可控精度算法介于上述两种算法之间,它能够使运算在可控的精度下进行计算.2.Matlab的数据显示格式,列表如下:表Matlab数据显示格式命令命令意义举例()formatshort短格式方式，显示5位定点十进制数3.1416formatlong长格式方式，显示15位定点十进制数formatshorte最优化短格式显示,5位加指数3.1416e+000formatlonge最优格式,15位加指数formatshortg5位定点或浮点格式3.1416formatlongg对双精度，显示15位定点或浮点格式，对单精度，显示7位定点或浮点格式formatshorteng至少5位加3位指数3.1416e+000formatlongeng16位加至少3位指数formathex十六进制格式方式400921fb54442d18formatbank银行格式.按元、角、分（小数点后具有两位）的固定格式3.14format++格式，以＋，—和空格分别表示中的正数，负数和零元素+format缺省时为默认短格式方式与formatshort相同3.1416formatrat分数格式形式.用有理数逼近显示数据355/113formatloose松散格式.数据之间有空行formatcompact紧凑格式.数据之间无空行vpa(date,n)将数据date以n位有效数字显示vpa(pi,5)=3.1416format并不影响matlab如何计算和存储变量的值.对浮点型变量的计算，即单精度或双精度，按合适的浮点精度进行，而不论变量是如何显示的.对整型变量采用整型数据.整型变量总是根据不同的类（class）以合适的数据位显示.3.Matlab的特殊变量ans:对最近输入的反应computer:当前计算机类型eps:浮点精度flops:计算浮点操作次数，现已不再常用i:虚部单位inf:无穷大inputname:输入参数名j:虚部单位nan:非数值nargin:输入参数的数目nargout:输出参数的数目(用户定义函数)pi:圆周率realmax:最大正浮点数realmin:最小正浮点数vararginvarargout:返回参数数目(matlab函数)cputime:CPU工作时间2线性方程组的直接解法2.1例题解答例2.1用Gauss消元法解方程组:解:直接建立求解该方程组的M文件Gauss.m如下:%求解例题2.1%高斯法求解线性方程组Ax=b%A为输入矩阵系数,b为方程组右端系数%方程组的解保存在x变量中%先输入方程系数A=[123;275;149];b=[16-3]';[m,n]=size(A);%检查系数正确性ifm~=nerror('矩阵A的行数和列数必须相同');return;endifm~=size(b)error('b的大小必须和A的行数或A的列数相同');return;end%再检查方程是否存在唯一解ifrank(A)~=rank([A,b])error('A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');return;end%这里采用增广矩阵行变换的方式求解c=n+1;A(:,c)=b;%%消元过程fork=1:n-1A(k+1:n,k:c)=A(k+1:n,k:c)-(A(k+1:n,k)/A(k,k))*A(k,k:c);End%%回代结果x=zeros(length(b),1);x(n)=A(n,c)/A(n,n);fork=n-1:-1:1x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k);end%显示计算结果disp('x=');disp(x);直接运行上面的M文件或在Matlab命令窗口中直接输入Gauss即可得出结果.在Matlab命令窗口中输入Gauss得出结果如下:>>Gaussx=2.00001.0000-1.0000扩展:Matlab求解线性方程的几种命令如下(方程组的一般形式可用矩阵和向量表示成:,但运用下列方法的前提必须保证所求解的方程为恰定方程,即方程组存在唯一的一组解):运用求逆思想:或;左除法,原理上是运用高斯消元法求解,但Matlab在实际执行过程中是通过分解法进行的(即先将矩阵A作分解,再回代计算):;用符号矩阵法计算,这种计算方法最接近精确值,但计算速度最慢:通过将矩阵施行初等行变换化成行简化阶梯形的办法,可以这样实现之:;.上面四种常用的办法示例如下:>>A=[123;275;149]%上面示例方程组系数A=123275149>>b=[16-3]'%方程组右端的系数b=16-3>>x1_1=inv(A)*b,x1_2=A^(-1)*b%方法一,求逆思想x1_1=2.00001.0000-1.0000x1_2=2.00001.0000-1.0000>>x2=A\b%方法二,左除思想x2=21-1>>x3=sym(A)\sym(b)%方法三,符号法x3=21-1>>C=[A,b],rref(C)%方法四,行简化阶梯形思想,最后输出结果的一列为解C=12312756149-3ans=10020101001-1例2.3用消元法思想求的逆阵.解:解法一:直接建立求解的M文件:Gauss_Jordan.m,源程序如下:%Gauss-Jordan法求例2.3clc;A=[1-10;223;-121];A1=A;%先保存原来的方阵A[n,m]=size(A);ifn~=merror('A必须为方阵');return;endA(:,n+1:2*n)=eye(n);%构造增广矩阵fork=1:n[l,m]=max(abs(A(k:n,k)));%按列选主元ifA(k+m-1,k)==0error('找到列最大的元素为零,错误');return;endifm~=1%交换Temp=A(k,:);A(k,:)=A(k+m-1,:);A(k+m-1,:)=Temp;endfori=1:nifi~=kA(i,:)=A(i,:)-A(k,:)*A(i,k)/A(k,k);endendendfori=n:(-1):1A(i,:)=A(i,:)/A(i,i);endA(:,1:n)=[];disp('A=');disp(A1);disp('用Gauss-Jandan算得矩阵A的逆矩阵为:');disp('inv(A)=');disp(A);clearTempiklmn;%清除临时变量在Matlab命令窗口中输入Gauss_Jordan回车后得到结果如下:A=1-10223-121用Gauss-Jandan算得矩阵A的逆矩阵为:inv(A)=-41-3-51-36-14解法二:用通过将矩阵施行初等行变换化成行简化阶梯形的办法,可以借助命令求解,非常方便:输入矩阵:>>A=[1-10;223;-121]A=1-10223-121>>C=[A,eye(length(A))]C=1-10100223010-121001>>invA=rref(C)invA=100-41-3010-51-30016-14后三列即为其逆阵,检验其正确性:>>c=invA(:,4:6)c=-41-3-51-36-14>>d=c*Ad=100010001可见求解正确.例2.4用分解LU的方法求解方程组解:解线性方程组中分解的L,U可以实现矩阵A的三角分解,使得A=L*U.L,U应该是下三角和上三角矩阵的,这样才利于回代求根.但是MATLAB中的LU分解与解线性方程组中的L,U不一样.MATLAB的LU分解命令调用格式为:[L,U]=lu(A)MATLAB计算出来的L是"准下三角"(交换L的行后才能成为真正的下三角阵),U为上三角矩阵,但它们还是满足A=L*U的.先录入矩阵系数.>>A=[2426;49615;26918;6151840]A=242649615269186151840>>b=[9232247]'b=9232247将A作分解,方法是使用矩阵分解的LU命令即可:>>[L,D]=lu(A)L=0.33331.0000-0.66671.00000.66671.0000000.3333-1.00001.000001.0000000U=6.000015.000018.000040.00000-1.0000-6.0000-11.666700-3.0000-7.0000000-0.3333再检验其正确性:>>C=L*UC=242649615269186151840解方程组>>y=L\by=47.0000-8.3333-2.00000.3333解方程组得到方程组的最终解:>>x=U\yx=0.50002.00003.0000-1.0000故方程组的最终解为:例2.5解方程组试用平方根法,改进的平方根法和追赶法分别解之.解:先输入矩阵:>>A=[610;141;0114]A=6101410114>>b=[624322]'b=624322平方根法:先对A矩阵作分解:>>L=chol(A)L=2.44950.4082001.95790.5108003.7066检验其正确性>>L'*Lans=6.00001.000001.00004.00001.000001.000014.0000将L转化为下三角矩阵:>>L=L'L=2.4495000.40821.9579000.51083.7066解方程组:>>y=L\by=2.449511.747385.2526再解方程组,得到最终解:>>x=L'\yx=1.0000-0.000023.0000故平方根法解上述方程的结果为:改进的平方根法:先对矩阵A作LDL分解:>>[L,D]=ldl(A)L=1.0000000.16671.0000000.26091.0000D=6.00000003.833300013.7391检验其分解正确性:>>L*D*L'ans=6101410114解方程组:>>y=L\by=623316解方程组:>>x=(D*L')\yx=1023故改进的平方根法解上述方程的结果亦为:追赶法:编制追赶法求解该方程的程序如下:%pursue.m%三对角线性方程组的追赶法解方程组例2.5%输入矩阵clc;A=[610;141;0114]f=[624322][n,m]=size(A);%分别取对角元素a=zeros(1,n);a(2:n)=diag(A,-1);c=diag(A,1);%此处用变量d存储A主对角线上的元素,因已用变量b存储方程右边的系数b=diag(A);ifb(1)==0error('主对角元素不能为0');return;end%初始计算,式(2.31)alpha(1)=b(1);beta(1)=c(1)/b(1);%按照公式(2.31)计算fori=2:n-1alpha(i)=b(i)-a(i)*beta(i-1);ifalpha(i)==0error('错误:在解方程过程中α为0');return;endbeta(i)=c(i)/alpha(i);end%对最后一行作计算alpha(n)=b(n)-a(n)*beta(n-1);ifalpha(n)==0error('错误:在解方程过程中最后一个α为0');return;end%以下按照公式(2.32)计算,解Ly=fy(1)=f(1)/b(1);fori=2:ny(i)=(f(i)-a(i)*y(i-1))/alpha(i);end%以下按照公式(2.33)计算,解Ux=yX(n)=y(n);fori=n-1:-1:1X(i)=y(i)-beta(i)*X(i+1);enddisp('X=');disp(X);在Matlab命令窗口输入pursue计算结果如下:>>pursueA=6101410114f=624322X=1023其中A为系数矩阵,f为矩阵右端的系数,最后计算结果为X.由以上计算可知追赶法解该方程的结果亦为:例2.6,求.解:输入矩阵:>>x=[10.5-0.3]x=1.00000.5000-0.3000计算其1-范数>>norm_1=norm(x,1)norm_1=1.8000计算其2-范数>>norm_2=norm(x)norm_2=1.1576计算其无穷大范数>>norm_inf=norm(x,inf)norm_inf=1还可以计算其无穷小范数(即各分量绝对值中的最小值)>>norm_minusInf=norm(x,-inf)norm_minusInf=0.3000例2.7求.解:先输入矩阵:>>A=[-210;1-21;01-2]A=-2101-2101-2求A的1-范数(列和范数):>>norm_1=norm(A,1)norm_1=4求解A的无穷大范数(行和范数):>>norm_inf=norm(A,inf)norm_inf=4求A的2-范数(最大特征值):>>norm_2=norm(A,2)norm_2=3.4142还可以求解A的F-范数:>>norm_F=norm(A,'fro')norm_F=4谱半径可以通过按其概念进行计算:对其特征值的绝对值取最大值即可:>>R_A=max(abs(eig(A)))R_A=3.4142例2.8n阶Hilbert矩阵考查其条件数.解:上述矩阵为抽象矩阵,而计算机只能进行有限次迭代,我们从考查其条件数,为此编制如下程序Hilb_cond10.m:%Hilb_cond10.m%考查从3阶至10的矩阵2—条件数fori=3:10h=hilb(i);condA(i)=cond(h,2);enddisp('ncond');fori=3:10s=sprintf('%d%f',i,condA(i));disp(s);end运行后得到如下结果:ncond3524.056778415513.7387395476607.250242614951058.6410057475367356.3703939493153986466.270940结果中左边为的阶数,右边为对应的条件数,从以上结果也可分析可知:当n的阶数稍高时,矩阵即出现严重的病态.例2.9求的条件数解:建立2阶矩阵:>>H=hilb(2)H=1.00000.50000.50000.3333求其2-条件数>>cond_2=cond(H,2)cond_2=19.2815求其1-条件数>>cond_1=cond(H,1)cond_1=27.0000求其无穷大条件数>>cond_inf=cond(H,inf)cond_inf=27.0000还可求其F-条件数>>cond_f=cond(H,'fro')cond_f=19.33332.2Matlab解线性方程组常用命令介绍1.求秩命令rank在解线性方程组时,为了判断是否存在解,我们应先判断矩阵的秩.调用格式为:c=rank(A)2.矩阵零空间向量null当线性方程组的秩时,方程组有无穷多个解,使用x=null(A)可得到满足的一个解向量,可为符号矩阵或数值矩阵:x=null(A)或x=null(sym(A))3.方阵的LU分解命令lu调用格式为:[L,U]=lu(A)L为准下三角矩阵,U为上三角矩阵,满足A=LU.[L,U,P]=lu(A)L为下三角矩阵,U为上三角矩阵,P为变换方阵元素位置的换位阵,满足PA=PL.其它调用格式请用helplu获得更多信息.4.Cholsky分解cholL=chol(A)其中L为一个下三角矩阵,满足.必须为正定矩阵.5.条件数condc=cond(A,p)A可为向量或矩阵,P取值为下列之一:1——向量或矩阵的返回1—条件数.2——返回向量或矩阵的2—条件数.inf——返回向量或矩阵的—条件数.'fro'——返回向量或矩阵的F—条件数.6.奇异值分解svd[U,S,V]=svd(A)将A分解为正交矩阵U,对角矩阵S和正交矩阵V的乘积,使得A=USVT.7.线性方程组的符号解linsolveX=linsolve(A,b)它等价于X=sym(A)\sym(b).返回方程组的符号解.3线性方程组的迭代解法3.1例题解答例3.1用Jacobi迭代法解以下方程:解:编制迭代计算的M文件程序如下:%Jacobi迭代法求解例3.1%A为方程组的增广矩阵clc;A=[10-2-13;-210-115;-1-2510]MAXTIME=50;%最多进行50次迭代eps=1e-5;%迭代误差[n,m]=size(A);x=zeros(n,1);%迭代初值y=zeros(n,1);k=0;%进入迭代计算disp('迭代过程X的值情况如下:')disp('X=');while1disp(x');fori=1:1:ns=0.0;forj=1:1:nifj~=is=s+A(i,j)*x(j);endy(i)=(A(i,n+1)-s)/A(i,i);endendfori=1:1:nmaxeps=max(0,abs(x(i)-y(i)));%检查是否满足迭代精度要求endifmaxeps<=eps%小于迭代精度退出迭代fori=1:1:nx(i)=y(i);%将结果赋给xendreturn;endfori=1:1:n%若不满足迭代精度要求继续进行迭代x(i)=y(i);y(i)=0.0;endk=k+1;ifk>MAXTIME%超过最大迭代次数退出error('超过最大迭代次数,退出');return;endend运行该程序结果如下:A=10-2-13-210-115-1-2510迭代过程X的值情况如下:X=0000.30001.50002.00000.80001.76002.66000.91801.92602.86400.97161.97002.95400.98941.98972.98230.99621.99612.99380.99861.99862.99770.99951.99952.99920.99981.99982.99970.99991.99992.99991.00002.00003.00001.00002.00003.0000容易看出迭代计算最后结果为:例3.2用Gauss-Seidel迭代法计算例3.1并作比较.解:编制求解程序Gauss_Seidel.m如下:%Gauss_Seidel.m%Gauss_Seidel迭代法求解例3.2%A为方程组的增广矩阵clc;formatlong;A=[10-2-13;-210-115;-1-2510][n,m]=size(A);%最多进行50次迭代Maxtime=50;%控制误差Eps=10E-5;%初始迭代值x=zeros(1,n);disp('x=');%迭代次数小于最大迭代次数,进入迭代fork=1:Maxtimedisp(x);fori=1:ns=0.0;forj=1:nifi~=js=s+A(i,j)*x(j);%计算和endendx(i)=(A(i,n+1)-s)/A(i,i);%求出此时迭代的值end%因为方程的精确解为整数,所以这里将迭代结果向整数靠近的误差作为判断迭代是否停止的条件ifsum((x-floor(x)).^2)<Epsbreak;end;end;X=x;disp('迭代结果:');Xformatshort;完成后直接在Matlab命令窗口中输入Gauss_Seidel回车后可得到如下结果:>>Gauss_SeidelA=10-2-13-210-115-1-2510x=0000.3000000000000001.5600000000000002.6840000000000000.8804000000000001.9444800000000002.9538720000000000.9842832000000001.9922438400000002.9937541760000000.9997021448448001.9998545228697602.9998822381168640.9999591283856381.9999800494888142.9999838454726540.9999943944450281.9999972634362712.9999977842635140.9999992311136061.9999996246490722.9999996960823500.9999998945380491.9999999485158452.9999999583139480.9999999855345641.9999999929383082.9999999942822360.9999999980158851.9999999990314012.9999999992157370.9999999997278541.9999999998671452.9999999998924290.9999999999626721.9999999999817782.9999999999852460.9999999999948801.9999999999975012.9999999999979760.9999999999992981.9999999999996572.9999999999997220.9999999999999041.9999999999999532.9999999999999620.9999999999999871.9999999999999942.9999999999999950.9999999999999981.9999999999999992.9999999999999991.0000000000000002.0000000000000003.000000000000000迭代结果:X=123可见对此题Gauss_Seidel法的收敛速度还是很快的.例3.3取,用Gauss-Seidel迭代法和SOR法计算下列方程组的解:解:Gauss-Seidel迭代法可利用上题中的程序,把输入矩阵A换掉就可以了,以下编制求解该程序的SOR.m%SOR法求解例3.3%w=1.45%方程组系数矩阵clc;A=[4-2-1;-24-2;-1-23]%方程组右端系数b=[0,-2,3]'w=1.45;%最大迭代次数Maxtime=100;%精度要求Eps=1E-5;%以15位小数显示formatlong;n=length(A);k=0;%初始迭代值x=ones(n,1);y=x;disp('迭代过程:');disp('x=');while1y=x;disp(x');%计算过程fori=1:ns=b(i);forj=1:nifj~=is=s-A(i,j)*x(j);endendifabs(A(i,i))<1E-10|k>=Maxtimeerror('已达最大迭代次数或矩阵系数近似为0,无法进行迭代');return;ends=s/A(i,i);x(i)=(1-w)*x(i)+w*s;endifnorm(y-x,inf)<Eps%达到精度要求退出计算break;endk=k+1;enddisp('最后迭代结果:');%最后的结果X=x'%设为默认显示格式formatshort;为了能有可比性,程序中的误差控制改用无穷大范数来度量(即将其改为:norm(y-x,inf)<Eps,并在循环中用y先保存迭代前x的值:y=x),此外,两个程序中的迭代初始值.以下为两种方法运行结果:(1)法运行结果如下:>>Gauss_SeidelA=4-2-10-24-2-2-1-233x=1110.7500000000000000.3750000000000001.5000000000000000.9322631785694630.9222812405563841.9256085532274100.9425427585850440.9340756559062271.9368980234658330.9586586469134330.9525664386408791.9545971747317300.9702