数值分析实验——何光辉

数值分析实验教程何光辉 重庆大学数学与统计学院二O一O年三月目 录0 Matlab介绍入门知识31 绪论161.1 例题解答161.2 Matlab中数值计算精度192 线性方程组的直接解法202.1 例题解答202.2 Matlab解线性方程组常用命令介绍343 线性方程组的迭代解法353.1 例题解答353.2 Matlab迭代解法用到的函数介绍484 方阵特征值和特征向量的计算484.1 例题解答484.2 Matlab关于方阵特征值为特征向量函数介绍565 非线性方程求根575.1 例题解答575.2 Matlab非线性方程求根的命令796 插值法796.1 例题解答796.2 Matlab插值函数介绍957 数据拟合和最佳平方逼近967.1 例题解答967.2 Matlab数据拟合命令介绍1068 数值积分与数值微分1078.1 例题解答1079 常微分方程数值解法1309.1 例题解答1309.2 Matlab常微分方程数值解常用命令介绍1460 Matlab介绍入门知识1. Matlab简介MATLAB的含义是矩阵实验室（MATRIX LABORATORY），主要用于方便矩阵的存取，其基本元素是无须定义维数的矩阵.MATLAB自问世以来,就是以数值计算称.MATLAB进行数值计算的基本单位是复数数组（或称阵列），这使得MATLAB高度“向量化”.经过十几年的完善和扩充，现已发展成为线性代数课程的标准工具.由于它不需定义数组的维数，并给出矩阵函数、特殊矩阵专门的库函数，使之在求解诸如信号处理、建模、系统识别、控制、优化等领域的问题时，显得大为简捷、高效、方便，这是其它高级语言所不能比拟的.MATLAB中包括了被称作工具箱（TOOLBOX）的各类应用问题的求解工具.工具箱实际上是对MATLAB进行扩展应用的一系列MATLAB函数（称为M文件），它可用来求解各类学科的问题，包括信号处理、图象处理、控制系统辨识、神经网络等.随着MATLAB版本的不断升级，其所含的工具箱的功能也越来越丰富，因此，应用范围也越来越广泛，MATLAB提供的工具箱已覆盖信号处理、系统控制、统计计算、优化计算、神经网络、小波分析、偏微分方程、模糊逻辑、动态系统模拟、系统辨识和符号运算等领域.当前它的使用范围涵盖了工业、电子、医疗、建筑等各行各业.MATLAB中包括了图形界面编辑GUI，让使用者也可以象VB、VC、VJ、DELPHI等那样进行一般的可视化的程序编辑.在命令窗口（matlab command window）键入simulink，就出现(SIMULINK) 窗口.以往十分困难的系统仿真问题，用SIMULINK只需拖动鼠标即可轻而易举地解决问题，这也是近来受到重视的原因所在. MATLAB 语言由美国 The MathWorks 开发,最早是由C.Moler用Fortran语言编写的，用来方便地调用LINPACK和EISPACK矩阵代数软件包的程序.后来他创立了MATHHWORKS公司，对MATLAB作了大量的、坚持不懈的改进.Cleve B.Moler是The MathWork公司的主席和首席科学家.曾任密歇系教授.他在两个计算机硬件制造商Intel公司的Hypercube组织和Arden Computers 公司工作了五年.他的主要专业兴趣在于数值分析和科学计算.他是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-asciiload 文件名 变量名表 -ascii其中，文件名可以带路径，但不需带扩展名.mat，命令隐含一定对.mat文件进行操作.变量名表中的变量个数不限，只要内存或文件中存在即可，变量名之间以空格分隔.当变量名表省略时，保存或装入全部变量.-ascii选项使文件以ASCII格式处理，省略该选项时文件将以二进制格式处理.save命令中的-append选项控制将变量追加到MAT文件中.(1) 向量的创建用步长生成法：数组=初值:步长(增量):终值 a=1:0.5:3a =1.0000 1.5000 2.0000 2.5000 3.0000用linspace生成：数组=linspace(初值,终值,等分点数目) b=linspace(1,3,5)b =1.0000 1.5000 2.0000 2.5000 3.0000列向量用分号（;）作为分行标记： c=1;2;3;4;c = 1 2 3 4若不想输出结果，在每一条语句后用分号作为结束符，若留空或用逗号结束，则在执行该语句后会把结果输出来. a+b; a+bans = 2 3 4 5 6(2) 矩阵的创建直接输入:最简单的建立矩阵的方法是从键盘直接输入矩阵的元素.具体方法如下：将矩阵的元素用方括号括起来，按矩阵行的顺序输入各元素，同一行的各元素之间用空格或逗号分隔，不同行的元素之间用分号分隔. A=1 2 3;4 5 6;2 3 5A = 1 2 3 4 5 6 2 3 5利用矩阵函数创建: B=magic(3)%魔方阵B = 8 1 6 3 5 7 4 9 2 C=hilb(3)%3阶Hilbert矩阵C = 1.0000 0.5000 0.3333 0.5000 0.3333 0.25000.3333 0.2500 0.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 = 9 3 9 7 10 13 6 12 7减法： B-Aans = 7 -1 3 -1 0 1 2 6 -3乘法： A*Bans = 26 38 26 71 83 7145 62 43除法: magic(3)/hilb(3)ans = 1.0e+003 * 0.2160 -1.1760 1.1400 0.0570 -0.4080 0.4500 -0.2280 1.2240 -1.1400在MATLAB中，有一种特殊的运算，因为其运算符是在有关算术运算符前面加点，所以叫点运算.点运算符有.*、./、.和.两矩阵进行点运算是指它们的对应元素进行相关运算，要求两矩阵的维参数相同. A.*Bans = 8 2 18 12 25 42 8 27 10MATLAB提供了6种关系运算符：(小于)、(大于)、=(大于或等于)、=(等于)、=(不等于). ABans = 0 1 0 1 0 0 0 0 1MATLAB提供了3种逻辑运算符：&(与)、|(或)和(非). 在逻辑运算中，确认非零元素为真，用1表示，零元素为假，用0表示.设参与逻辑运算的是两个标量a和b，那么， a&b a,b全为非零时，运算结果为1，否则为0. a|b a,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矩阵第ii+m行的全部元素；A(:,k:k+m)表示取A矩阵第kk+m列的全部元素，A(i:i+m,k:k+m)表示取A矩阵第ii+m行内，并在第kk+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 = 8 1 6 3 5 7 4 9 2 B=inv(A)B = 0.1472 -0.1444 0.0639 -0.0611 0.0222 0.1056 -0.0194 0.1889 -0.1028 A*Bans = 1.0000 0 -0.0000 -0.0000 1.0000 00.0000 0 1.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)绘制出具有不同纵坐标标度的两个图形.hold on/off保持原有图形还是刷新原有图形,不带参数的hold命令在两种状态之间进行切换.plot(x1,y1,选项1,x2,y2,选项2,xn,yn,选项n)设置曲线样式进行绘图.选项字段见下表:表 Matlab常用线形与颜色标记表符号线型符号线型符号线型颜色含义-实线.实点标记朝上三角g绿色-虚线o圆圈朝右三角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) 坐标控制axis axis(xmin xmax ymin ymax zmin zmax)axis函数功能丰富，常用的格式还有：axis equal：纵、横坐标轴采用等长刻度.axis square：产生正方形坐标系(缺省为矩形).axis auto：使用缺省设置.axis off：取消坐标轴.axis on：显示坐标轴.grid on/off:网格开/关box on/off:加/不加边框线上述命令示例如下: x=1:length(peaks); plot(x,peaks); box on; 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*-1 1 -1 1);(5) 二维符号函数曲线专用命令ezplotf = f(x)时:ezplot(f)：在默认区间-2x2绘制f = f(x)的图形.ezplot(f, a,b)：在区间axb绘制f = f(x)的图形f = f(x,y)时:ezplot(f)：在默认区间-2x2和-2y2绘制f(x,y) = 0的图形.ezplot(f, xmin,xmax,ymin,ymax)：在区间xminxxmax和yminyymax绘制f(x,y) = 0的图形ezplot(f, a,b)：在区间axb和ay b绘制f(x,y) = 0的图形若x = x(t),y = y(t):ezplot(x,y)：在默认区间0t2绘制x=x(t)和y=y(t)的图形.ezplot(x,y, tmin,tmax)：在区间tmin t figure;ezplot(cos(tan(pi*x), 0,1);(6) 图形窗口的分割subplotsubplot(m,n,p)该函数将当前图形窗口分成mn个绘图区，即每行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 条件 语句组1 else 语句组2 end当条件成立时，执行语句组1，否则执行语句组2，语句组1或语句组2执行后，再执行if语句的后继语句.(6) 多分支if语句：if 条件1 语句组1 elseif 条件2 语句组2 elseif 条件m 语句组m else 语句组n end语句用于实现多分支选择结构.(7)switch语句: switch 表达式 case 表达式1 语句组1 case 表达式2 语句组2 case 表达式m 语句组m otherwise 语句组n end(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 计算,.解:创建符号函数: syms x; f=sym(sin(x) f = sin(x)展开至7阶泰勒级数： h=taylor(f,8,0) h = x-1/6*x3+1/120*x5-1/5040*x7求泰勒级数在处的函数值： subs(h,x,0.5) ans = 0.4127也可以通过内联函数来求解:H=inline(h) H = Inline function: H(x) = x-1./6.*x.3+1./120.*x.5-1./5040.*x.7 feval(H,0.5) ans = 0.4127例 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;for i=1:5 I=I+(f(i)+f(i+1)/2*0.2; %将每一小曲边的梯形累加起来作为积分值End vpa(I,9) %取结果的小数精度为9位小数 ans = .例 1.3略例 1.4 不用开平方根计算的值.解:解法一(符号法): A=sym(a); sqrt(A) ans = a(1/2)解法二(数值法):按以下迭代公式迭代计算近似值:建立函数文件msqrt.mfunction x=msqrt(x0,a)%用迭代法近似计算平方根%x0为初始迭代值,a为开平方数 format long;x=zeros(20,1);x(1)=x0;for i=2:20 x(i)=1/2*(x(i-1)+a/x(i-1);enddisp(x);用编写的函数计算,: msqrt(2,3); 2.0000 1.0000 1.2857 1.4727 1.8877 1.8877 1.8877 1.8877 1.8877 1.8877 1.8877 1.8877 1.8877 1.8877 1.8877 1.8877 1.8877 1.8877 1.8877 1.8877上述结果为迭代过程计算的中间结果,分析数据可知迭代收敛速度快,只需四次计算即可计算出较为准确的数值.例 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-b a = 1/759 b = 1/760 c = 1/可见符号运算准确但耗费运算时间.例 1.7 略.例 1.8 解方程.解:符号法解方程: x=solve(x2-18*x+1,x) x = 9+4*5(1/2) 9-4*5(1/2)将结果保留小数点6位: vpa(x,6) ans = 17.9443 .5572e-11.2 Matlab中数值计算精度1. Matlab中有三种运算精度,它们分别为数值算法、符号算法和可控精度算法,将它们分别介绍如下:(1) 数值算法把每个数取为16位,计算按浮点运算进行,它是运算速度最快的一种算法.(2) 符号算法把每个数都变为符号量,运算按有理量计算进行,它的优点是能够得到精确结果,缺点是占用空间大,并且运算速度最慢.(3) 可控精度算法介于上述两种算法之间,它能够使运算在可控的精度下进行计算.2. Matlab的数据显示格式,列表如下: 表 Matlab数据显示格式命令命令意义举例()format short短格式方式，显示5位定点十进制数3.1416format long长格式方式，显示15位定点十进制数3.9793format short e最优化短格式显示,5位加指数3.1416e+000format long e最优格式,15位加指数3.9793e+000format short g5位定点或浮点格式3.1416format long g对双精度，显示15位定点或浮点格式，对单精度，显示7位定点或浮点格式3.979format short eng至少5位加3位指数3.1416e+000format long eng16位加至少3位指数3.979e+000format hex十六进制格式方式fb54442d18format bank银行格式.按元、角、分（小数点后具有两位）的固定格式3.14format +格式，以，和空格分别表示中的正数，负数和零元素+format缺省时为默认短格式方式与format short相同3.1416format rat分数格式形式.用有理数逼近显示数据355/113 format loose松散格式.数据之间有空行format compact紧凑格式.数据之间无空行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 :最小正浮点数 varargin varargout :返回参数数目(matlab函数) cputime:CPU工作时间2 线性方程组的直接解法2.1 例题解答例 2.1 用Gauss消元法解方程组:解:直接建立求解该方程组的M文件Gauss.m如下:%求解例题2.1%高斯法求解线性方程组Ax=b%A为输入矩阵系数,b为方程组右端系数%方程组的解保存在x变量中%先输入方程系数A=1 2 3;2 7 5;1 4 9;b=1 6 -3;m,n=size(A);%检查系数正确性if m=n error(矩阵A的行数和列数必须相同); return;endif m=size(b) error(b的大小必须和A的行数或A的列数相同); return;end%再检查方程是否存在唯一解if rank(A)=rank(A,b) error(A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解); return;end%这里采用增广矩阵行变换的方式求解c=n+1;A(:,c)=b; %消元过程for k=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);for k=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.0000 1.0000 -1.0000扩展:Matlab求解线性方程的几种命令如下(方程组的一般形式可用矩阵和向量表示成: ,但运用下列方法的前提必须保证所求解的方程为恰定方程,即方程组存在唯一的一组解):运用求逆思想: 或 ;左除法,原理上是运用高斯消元法求解,但Matlab在实际执行过程中是通过分解法进行的(即先将矩阵A作分解,再回代计算): ;用符号矩阵法计算,这种计算方法最接近精确值,但计算速度最慢: 通过将矩阵施行初等行变换化成行简化阶梯形的办法,可以这样实现之: ; .上面四种常用的办法示例如下: A=1 2 3;2 7 5;1 4 9 %上面示例方程组系数A = 1 2 3 2 7 5 1 4 9 b=1 6 -3 %方程组右端的系数b = 1 6-3 x1_1=inv(A)*b,x1_2=A(-1)*b %方法一,求逆思想x1_1 = 2.0000 1.0000 -1.0000x1_2 = 2.0000 1.0000 -1.0000 x2=Ab %方法二,左除思想x2 = 2 1-1 x3=sym(A)sym(b) %方法三,符号法 x3 = 2 1 -1 C=A,b,rref(C) %方法四, 行简化阶梯形思想,最后输出结果的一列为解C = 1 2 3 1 2 7 5 6 1 4 9 -3ans = 1 0 0 2 0 1 0 1 0 0 1 -1例 2.2 用选列主元的Gauss消元法解如下方程:解:直接建立求解该方程的M文件Gauss_line.m,求解程序编制如下:%求解例题2.2%高斯列主元消元法求解线性方程组Ax=b%A为输入矩阵系数,b为方程组右端系数%方程组的解保存在x变量中format long;%设置为长格式显示,显示15位小数A=0.00001,1;2,1b=1.00001,3m,n=size(A);%先检查系数正确性if m=n error(矩阵A的行数和列数必须相同); return;endif m=size(b) error(b的大小必须和A的行数或A的列数相同); return;end%再检查方程是否存在唯一解if rank(A)=rank(A,b) error(A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解); return;endc=n+1;A(:,c)=b; %(增广)for k=1:n-1r,m=max(abs(A(k:n,k); %选主元 m=m+k-1; %修正操作行的值 if(A(m,k)=0) if(m=k) A(k m,:)=A(m k,:); %换行 end A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k)*A(k, k:c); %消去 endendx=zeros(length(b),1); %回代求解x(n)=A(n,c)/A(n,n);for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n)/A(k,k);enddisp(X=);disp(x);format short;%设置为默认格式显示,显示5位运行得到输出结果如下:A = 0.0000 1.0000 2.0000 1.0000b = 1.0000 3.0000X= 1.0000 1.0000例 2.3 用消元法思想求的逆阵.解:解法一:直接建立求解的M文件:Gauss_Jordan.m,源程序如下:%Gauss-Jordan法求例2.3clc;A=1 -1 0;2 2 3;-1 2 1;A1=A;%先保存原来的方阵An,m=size(A); if n=m error(A必须为方阵); return; endA(:,n+1:2*n)=eye(n);%构造增广矩阵for k=1:n l,m=max(abs(A(k:n,k);%按列选主元 if A(k+m-1,k)=0 error(找到列最大的元素为零,错误); return; end if m=1 %交换 Temp=A(k,:); A(k,:)=A(k+m-1,:); A(k+m-1,:)=Temp; end for i=1:n if i=k A(i,:)=A(i,:)-A(k,:)*A(i,k)/A(k,k); end end endfor i=n:(-1):1 A(i,:)=A(i,:)/A(i,i);endA(:,1:n)=;disp(A=);disp(A1);disp(用Gauss-Jandan算得矩阵A的逆矩阵为:);disp(inv(A)=);disp(A); clear Temp i k l m n;%清除临时变量在Matlab命令窗口中输入Gauss_Jordan回车后得到结果如下:A= 1 -1 0 2 2 3 -1 2 1用Gauss-Jandan算得矩阵A的逆矩阵为:inv(A)= -4 1 -3 -5 1 -3 6 -1 4解法二:用通过将矩阵施行初等行变换化成行简化阶梯形的办法,可以借助命令求解,非常方便:输入矩阵: A=1 -1 0;2 2 3;-1 2 1A = 1 -1 0 2 2 3-1 2 1 C=A,eye(length(A)C = 1 -1 0 1 0 0 2 2 3 0 1 0-1 2 1 0 0 1 invA=rref(C)invA = 1 0 0 -4 1 -3 0 1 0 -5 1 -3 0 0 1 6 -1 4后三列即为其逆阵,检验其正确性: c=invA(:,4:6)c = -4 1 -3 -5 1 -3 6 -1 4 d=c*Ad = 1 0 0 0 1 0 0 0 1可见求解正确.例 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=2 4 2 6;4 9 6 15;2 6 9 18;6 15 18 40A = 2 4 2 6 4 9 6 15 2 6 9 18 6 15 18 40 b=9 23 22 47b = 9 23 22 47将A作分解,方法是使用矩阵分解的LU命令即可:L,D=lu(A)L = 0.3333 1.0000 -0.6667 1.0000 0.6667 1.0000 0 0 0.3333 -1.0000 1.0000 0 1.0000 0 0 0U = 6.0000 15.0000 18.0000 40.0000 0 -1.0000 -6.0000 -11.6667 0 0 -3.0000 -7.0000 0 0 0 -0.3333再检验其正确性:C=L*U C = 2 4 2 6 4 9 6 15 2 6 9 18 6 15 18 40解方程组y=Lb y = 47.0000 -8.3333 -2.0000 0.3333 解方程组得到方程组的最终解:x=Uy x = 0.5000 2.0000 3.0000 -1.0000 故方程组的最终解为: 例 2.5 解方程组试用平方根法,改进的平方根法和追赶法分别解之.解:先输入矩阵:A=6 1 0;1 4 1;0 1 14 A = 6 1 0 1 4 1 0 1 14 b=6 24 322 b = 6 24