计算方法数值分析实验

数值分析实验董海云数理学院数学实验教学中心2008年3月目录0MATLAB介绍入门知识31绪论1611例题解答1612MATLAB中数值计算精度192线性方程组的直接解法2021例题解答2022MATLAB解线性方程组常用命令介绍343线性方程组的迭代解法3531例题解答3532MATLAB迭代解法用到的函数介绍484方阵特征值和特征向量的计算4841例题解答4842MATLAB关于方阵特征值为特征向量函数介绍565非线性方程求根5751例题解答5752MATLAB非线性方程求根的命令796插值法7961例题解答7962MATLAB插值函数介绍957数据拟合和最佳平方逼近9671例题解答9672MATLAB数据拟合命令介绍1068数值积分与数值微分10781例题解答1079常微分方程数值解法13091例题解答13092MATLAB常微分方程数值解常用命令介绍1460MATLAB介绍入门知识1MATLAB简介MATLAB的含义是矩阵实验室（MATRIXLABORATORY），主要用于方便矩阵的存取，其基本元素是无须定义维数的矩阵MATLAB自问世以来,就是以数值计算称MATLAB进行数值计算的基本单位是复数数组（或称阵列），这使得MATLAB高度“向量化”经过十几年的完善和扩充，现已发展成为线性代数课程的标准工具由于它不需定义数组的维数，并给出矩阵函数、特殊矩阵专门的库函数，使之在求解诸如信号处理、建模、系统识别、控制、优化等领域的问题时，显得大为简捷、高效、方便，这是其它高级语言所不能比拟的MATLAB中包括了被称作工具箱（TOOLBOX）的各类应用问题的求解工具工具箱实际上是对MATLAB进行扩展应用的一系列MATLAB函数（称为M文件），它可用来求解各类学科的问题，包括信号处理、图象处理、控制系统辨识、神经网络等随着MATLAB版本的不断升级，其所含的工具箱的功能也越来越丰富，因此，应用范围也越来越广泛，MATLAB提供的工具箱已覆盖信号处理、系统控制、统计计算、优化计算、神经网络、小波分析、偏微分方程、模糊逻辑、动态系统模拟、系统辨识和符号运算等领域当前它的使用范围涵盖了工业、电子、医疗、建筑等各行各业MATLAB中包括了图形界面编辑GUI，让使用者也可以象VB、VC、VJ、DELPHI等那样进行一般的可视化的程序编辑在命令窗口（MATLABCOMMANDWINDOW）键入SIMULINK，就出现SIMULINK窗口以往十分困难的系统仿真问题，用SIMULINK只需拖动鼠标即可轻而易举地解决问题，这也是近来受到重视的原因所在MATLAB语言由美国THEMATHWORKS开发,最早是由CMOLER用FORTRAN语言编写的，用来方便地调用LINPACK和EISPACK矩阵代数软件包的程序后来他创立了MATHHWORKS公司，对MATLAB作了大量的、坚持不懈的改进CLEVEBMOLER是THEMATHWORK公司的主席和首席科学家曾任密歇系教授他在两个计算机硬件制造商INTEL公司的HYPERCUBE组织和ARDENCOMPUTERS公司工作了五年他的主要专业兴趣在于数值分析和科学计算他是MATLAB软件的创始者，也是著名的矩阵计算软件包LINPACK和EISPACK的著作这一，已撰写了三本有相关数值方法的教材同时，他在SIAM（美国工业与应用数学学会）历任期刊编辑、委员会成员和副总裁，并从1996年开始担任理事会成员2MATLAB入门知识MATLAB变量名是以字母开头，后接字母、数字或下划线的字符序列，最多63个字符在MATLAB中，变量名区分字母的大小写赋值语句变量表达式或表达式其中表达式是用运算符将有关运算量连接起来的式子，其结果是一个矩阵CLEAR命令用于删除MATLAB工作空间中的变量WHO和WHOS这两个命令用于显示在MATLAB工作空间中已经驻留的变量名清单WHO命令只显示出驻留变量的名称，WHOS在给出变量名的同时，还给出它们的大小、所占字节数及数据类型等信息利用MAT文件可以把当前MATLAB工作空间中的一些有用变量长久地保留下来，扩展名是MATMAT文件的生成和装入由SAVE和LOAD命令来完成常用格式为SAVE文件名变量名表APPENDASCIILOAD文件名变量名表ASCII其中，文件名可以带路径，但不需带扩展名MAT，命令隐含一定对MAT文件进行操作变量名表中的变量个数不限，只要内存或文件中存在即可，变量名之间以空格分隔当变量名表省略时，保存或装入全部变量ASCII选项使文件以ASCII格式处理，省略该选项时文件将以二进制格式处理SAVE命令中的APPEND选项控制将变量追加到MAT文件中1向量的创建用步长生成法数组初值步长增量终值A1053A1000015000200002500030000用LINSPACE生成数组LINSPACE初值,终值,等分点数目BLINSPACE1,3,5B1000015000200002500030000列向量用分号（）作为分行标记C1234C1234若不想输出结果，在每一条语句后用分号作为结束符，若留空或用逗号结束，则在执行该语句后会把结果输出来ABABANS234562矩阵的创建直接输入最简单的建立矩阵的方法是从键盘直接输入矩阵的元素具体方法如下将矩阵的元素用方括号括起来，按矩阵行的顺序输入各元素，同一行的各元素之间用空格或逗号分隔，不同行的元素之间用分号分隔A123456235A123456235利用矩阵函数创建BMAGIC3魔方阵B816357492CHILB33阶HILBERT矩阵C100000500003333050000333302500033330250002000MATLAB中用引导注释其它创建矩阵函数还有EYEM,N生成M行N列单位矩阵ZEROSM,N生成M行N列全零矩阵ONESM,N生成全1矩阵RANDM,N生成随机矩阵RAND生成一个随机数DIAGA取A的对角线元素TRILA取A的下三角元素TRIUA取A的上三角元素HILBN生成N维HILBERT矩阵RANDNN产生均值为0，方差为1的标准正态分布随机矩阵VANDERV生成以向量V为基础向量的范得蒙矩阵INVHILBN求N阶的希尔伯特矩阵的逆矩阵TOEPLITZX,Y生成一个以X为第一列，Y为第一行的托普利兹矩阵COMPANP生成伴随矩阵,P是一个多项式的系数向量，高次幂系数排在前，低次幂排在后PASCALN生成一个N阶帕斯卡矩阵COMPAN生成伴随矩阵3矩阵运算MATLAB的基本算术运算有加、减、乘、/右除、左除、乘方加法ABANS939710136127减法BAANS713101263乘法ABANS263826718371456243除法MAGIC3/HILB3ANS10E003021601176011400005700408004500022801224011400在MATLAB中，有一种特殊的运算，因为其运算符是在有关算术运算符前面加点，所以叫点运算点运算符有、/、和两矩阵进行点运算是指它们的对应元素进行相关运算，要求两矩阵的维参数相同ABANS821812254282710MATLAB提供了6种关系运算符大于、大于或等于、等于、不等于ABANS010100001MATLAB提供了3种逻辑运算符PLOTX,PEAKSBOXONTITLE绘制混合图形XLABELX轴YLABELY轴绘制图像为0510152025303540455086420246810一一一一一一X一Y一4二维数值函数的专用绘图函数FPLOTFPLOTFUNCTIONNAME,A,B,TOL,选项其中FUNCTIONNAME为函数名，以字符串形式出现，A,B为绘图区间，TOL为相对允许误差，其系统默认值为2E3选项定义与PLOT函数相同FPLOTXTANX,SINX,COSX,2PI1111642024664202465二维符号函数曲线专用命令EZPLOTFFX时EZPLOTF在默认区间2FIGUREEZPLOTCOSTANPIX,0,1001020304050607080911050051XCOSTANX6图形窗口的分割SUBPLOTSUBPLOTM,N,P该函数将当前图形窗口分成MN个绘图区，即每行N个，共M行，区号按行优先编号，且选定第P个区为当前活动区在每一个绘图区允许以不同的坐标系单独绘制图形7其他坐标系下的二维数据曲线图对数坐标图形SEMILOGXX1,Y1,选项1,X2,Y2,选项2,SEMILOGYX1,Y1,选项1,X2,Y2,选项2,LOGLOGX1,Y1,选项1,X2,Y2,选项2,极坐标图POLARPOLARTHETA,R,选项其中THETA为极坐标极角，R为极坐标矢径，选项的内容与PLOT函数相似二维统计分析图BARX,Y,选项条形图STAIRSX,Y,选项阶梯图STEMX,Y,选项杆图FILLX1,Y1,选项1,X2,Y2,选项2,填充图8三维曲线PLOT3PLOT3X1,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对应列元素绘制三维曲线，曲线条数等于矩阵列数T0018PIPLOT3SINT,COST,T105005110500510510152025309产生三维数据在MATLAB中，利用MESHGRID函数产生平面区域内的网格坐标矩阵其格式为X,YMESHGRIDX,Y语句执行后，矩阵X的每一行都是向量X，行数等于向量Y的元素的个数，矩阵Y的每一列都是向量Y，列数等于向量X的元素的个数10绘制三维曲面的函数SURF函数和MESH函数的调用格式为MESHX,Y,Z,CSURFX,Y,Z,C一般情况下，X,Y,Z是维数相同的矩阵X,Y是网格坐标矩阵，Z是网格点上的高度矩阵，C用于指定在不同高度下的颜色范围11标准三维曲面SPHERE函数的调用格式为X,Y,ZSPHERENCYLINDER函数的调用格式为X,Y,ZCYLINDERR,NMATLAB还有一个PEAKS函数，称为多峰函数，常用于三维曲面的演示12其他三维绘图指令介绍BAR3函数绘制三维条形图，常用格式为BAR3YBAR3X,YSTEM3函数绘制离散序列数据的三维杆图，常用格式为STEM3ZSTEM3X,Y,ZPIE3函数绘制三维饼图，常用格式为PIE3XFILL3函数等效于三维函数FILL，可在三维空间内绘制出填充过的多边形，常用格式为FILL3X,Y,Z,C5程序控制结构1数据的输入AINPUT提示信息，选项其中提示信息为一个字符串，用于提示用户输入什么样的数据如果在INPUT函数调用时采用S选项，则允许用户输入一个字符串2数据的输出DISP输出项3程序的暂停PAUSE延迟秒数如果省略延迟时间，直接使用PAUSE，则将暂停程序，直到用户按任一键后程序继续执行若要强行中止程序的运行可使用CTRLC命令4单分支IF语句IF条件语句组END当条件成立时，则执行语句组，执行完之后继续执行IF语句的后继语句，若条件不成立，则直接执行IF语句的后继语句5双分支IF语句IF条件语句组1ELSE语句组2END当条件成立时，执行语句组1，否则执行语句组2，语句组1或语句组2执行后，再执行IF语句的后继语句6多分支IF语句IF条件1语句组1ELSEIF条件2语句组2ELSEIF条件M语句组MELSE语句组NEND语句用于实现多分支选择结构7SWITCH语句SWITCH表达式CASE表达式1语句组1CASE表达式2语句组2CASE表达式M语句组MOTHERWISE语句组NEND8TRY语句语句格式为TRY语句组1CATCH语句组2ENDTRY语句先试探性执行语句组1，如果语句组1在执行过程中出现错误，则将错误信息赋给保留的LASTERR变量，并转去执行语句组29FOR语句FOR语句的格式为FOR循环变量表达式1表达式2表达式3循环体语句END其中表达式1的值为循环变量的初值，表达式2的值为步长，表达式3的值为循环变量的终值步长为1时，表达式2可以省略FOR语句更一般的格式为FOR循环变量矩阵表达式循环体语句END执行过程是依次将矩阵的各列元素赋给循环变量，然后执行循环体语句，直至各列元素处理完毕10WHILE语句WHILE语句的一般格式为WHILE条件循环体语句END其执行过程为若条件成立，则执行循环体语句，执行后再判断条件是否成立，如果不成立则跳出循环11BREAK语句和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绪论11例题解答例11计算,SINX04解创建符号函数SYMSXFSYMSINXFSINX展开至7阶泰勒级数HTAYLORF,8,0HX1/6X31/120X51/5040X7求泰勒级数在处的函数值05XSUBSH,X,05ANS0479425533234127也可以通过内联函数来求解HINLINEHHINLINEFUNCTIONHXX1/6X31/120X51/5040X7FEVALH,05ANS0479425533234127例12计算积分值10IDX解解法一符号法IINT1/1X,X,0,1ILOG2解法二数值法X0021将0,1等分为4等份F1/1X分别计算每一个等分点的函数值I0FORI15IIFIFI1/202将每一小曲边的梯形累加起来作为积分值ENDVPAI,9取结果的小数精度为9位小数ANS695634921例13略例14不用开平方根计算的值0A解解法一符号法ASYMASQRTAANSA1/2解法二数值法按以下迭代公式迭代计算近似值1,01,22KKAXX建立函数文件MSQRTMFUNCTIONXMSQRTX0,A用迭代法近似计算平方根X0为初始迭代值,A为开平方数FORMATLONGXZEROS20,1X1X0FORI220XI1/2XI1A/XI1ENDDISPX用编写的函数计算,302XMSQRT2,320000000000000001750000000000000173214285714285717320508100147271732050807568877173205080756887717320508075688771732050807568877173205080756887717320508075688771732050807568877173205080756887717320508075688771732050807568877173205080756887717320508075688771732050807568877173205080756887717320508075688771732050807568877上述结果为迭代过程计算的中间结果,分析数据可知迭代收敛速度快,只需四次计算即可计算出较为准确的数值例15略例16计算,视已知数为精确数,用4位浮点数计算175960解直接在MATLAB中输入式子1/7591/760ANS17336E006若先转化为浮点数再运算可得A1/759,B1/760,ABA00013B00013ANS17336E006可见MATLBA在计算时,数据结构都取为双精度而提高了运算准确度若以符号运算计算之,有ASYM1/759,BSYM1/760,CABA1/759B1/760C1/576840可见符号运算准确但耗费运算时间例17略例18解方程2180X解符号法解方程XSOLVEX218X1,XX9451/29451/2将结果保留小数点6位VPAX,6ANS1794435572E112MATLAB中数值计算精度1MATLAB中有三种运算精度,它们分别为数值算法、符号算法和可控精度算法,将它们分别介绍如下1数值算法把每个数取为16位,计算按浮点运算进行,它是运算速度最快的一种算法2符号算法把每个数都变为符号量,运算按有理量计算进行,它的优点是能够得到精确结果,缺点是占用空间大,并且运算速度最慢3可控精度算法介于上述两种算法之间,它能够使运算在可控的精度下进行计算2MATLAB的数据显示格式,列表如下表MATLAB数据显示格式命令命令意义举例FORMATSHORT短格式方式，显示5位定点十进制数31416FORMATLONG长格式方式，显示15位定点十进制数3141592653589793FORMATSHORTE最优化短格式显示,5位加指数31416E000FORMATLONGE最优格式,15位加指数3141592653589793E000FORMATSHORTG5位定点或浮点格式31416FORMATLONGG对双精度，显示15位定点或浮点格式，对单精度，显示7位定点或浮点格式314159265358979FORMATSHORTENG至少5位加3位指数31416E000FORMATLONGENG16位加至少3位指数314159265358979E000FORMATHEX十六进制格式方式400921FB54442D18FORMATBANK银行格式按元、角、分（小数点后具有两位）的固定格式314FORMAT格式，以，和空格分别表示中的正数，负数和零元素FORMAT缺省时为默认短格式方式与31416FORMATSHORT相同FORMATRAT分数格式形式用有理数逼近显示数据355/113FORMATLOOSE松散格式数据之间有空行FORMATCOMPACT紧凑格式数据之间无空行VPADATE,N将数据DATE以N位有效数字显示VPAPI,531416FORMAT并不影响MATLAB如何计算和存储变量的值对浮点型变量的计算，即单精度或双精度，按合适的浮点精度进行，而不论变量是如何显示的对整型变量采用整型数据整型变量总是根据不同的类（CLASS）以合适的数据位显示3MATLAB的特殊变量ANS对最近输入的反应COMPUTER当前计算机类型EPS浮点精度FLOPS计算浮点操作次数，现已不再常用I虚部单位INF无穷大INPUTNAME输入参数名J虚部单位NAN非数值NARGIN输入参数的数目NARGOUT输出参数的数目用户定义函数PI圆周率REALMAX最大正浮点数REALMIN最小正浮点数VARARGINVARARGOUT返回参数数目MATLAB函数CPUTIMECPU工作时间2线性方程组的直接解法21例题解答例21用GAUSS消元法解方程组123175649XX解直接建立求解该方程组的M文件GAUSSM如下求解例题21高斯法求解线性方程组AXBA为输入矩阵系数,B为方程组右端系数方程组的解保存在X变量中先输入方程系数A123275149B163M,NSIZEA检查系数正确性IFMNERROR矩阵A的行数和列数必须相同RETURNENDIFMSIZEBERRORB的大小必须和A的行数或A的列数相同RETURNEND再检查方程是否存在唯一解IFRANKARANKA,BERRORA矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解RETURNEND这里采用增广矩阵行变换的方式求解CN1A,CB消元过程FORK1N1AK1N,KCAK1N,KCAK1N,K/AK,KAK,KCEND回代结果XZEROSLENGTHB,1XNAN,C/AN,NFORKN111XKAK,CAK,K1NXK1N/AK,KEND显示计算结果DISPXDISPX直接运行上面的M文件或在MATLAB命令窗口中直接输入GAUSS即可得出结果在MATLAB命令窗口中输入GAUSS得出结果如下GAUSSX200001000010000扩展MATLAB求解线性方程的几种命令如下方程组的一般形式可用矩阵和向量表示成,但运用下列方法的前提必须保证所求解的方程为恰定方程,即方程AXB组存在唯一的一组解运用求逆思想或1XINVAB1XB左除法,原理上是运用高斯消元法求解,但MATLAB在实际执行过程中是通过分解法进行的即先将矩阵A作分解,再回代计算LULUXA用符号矩阵法计算,这种计算方法最接近精确值,但计算速度最慢XSYMAB通过将矩阵施行初等行变换化成行简化阶梯形的办法,可以这样实现之,CREF上面四种常用的办法示例如下A123275149上面示例方程组系数A123275149B163方程组右端的系数B163X1_1INVAB,X1_2A1B方法一,求逆思想X1_1200001000010000X1_2200001000010000X2AB方法二,左除思想X2211X3SYMASYMB方法三,符号法X3211CA,B,RREFC方法四,行简化阶梯形思想,最后输出结果的一列为解C123127561493ANS100201010011例22用选列主元的GAUSS消元法解如下方程120013X解直接建立求解该方程的M文件GAUSS_LINEM,求解程序编制如下求解例题22高斯列主元消元法求解线性方程组AXBA为输入矩阵系数,B为方程组右端系数方程组的解保存在X变量中FORMATLONG设置为长格式显示,显示15位小数A000001,12,1B100001,3M,NSIZEA先检查系数正确性IFMNERROR矩阵A的行数和列数必须相同RETURNENDIFMSIZEBERRORB的大小必须和A的行数或A的列数相同RETURNEND再检查方程是否存在唯一解IFRANKARANKA,BERRORA矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解RETURNENDCN1A,CB增广FORK1N1R,MMAXABSAKN,K选主元MMK1修正操作行的值IFAM,K0IFMKAKM,AMK,换行ENDAK1N,KCAK1N,KCAK1N,K/AK,KAK,KC消去ENDENDXZEROSLENGTHB,1回代求解XNAN,C/AN,NFORKN111XKAK,CAK,K1NXK1N/AK,KENDDISPXDISPXFORMATSHORT设置为默认格式显示,显示5位运行得到输出结果如下A0000010000000000100000000000000020000000000000001000000000000000B10000100000000003000000000000000X10000000000000001000000000000000例23用消元法思想求GAUSJORDN1023A的逆阵解解法一直接建立求解的M文件GAUSS_JORDANM,源程序如下GAUSSJORDAN法求例23CLCA110223121A1A先保存原来的方阵AN,MSIZEAIFNMERRORA必须为方阵RETURNENDA,N12NEYEN构造增广矩阵FORK1NL,MMAXABSAKN,K按列选主元IFAKM1,K0ERROR找到列最大的元素为零,错误RETURNENDIFM1交换TEMPAK,AK,AKM1,AKM1,TEMPENDFORI1NIFIKAI,AI,AK,AI,K/AK,KENDENDENDFORIN11AI,AI,/AI,IENDA,1NDISPADISPA1DISP用GAUSSJANDAN算得矩阵A的逆矩阵为DISPINVADISPACLEARTEMPIKLMN清除临时变量在MATLAB命令窗口中输入GAUSS_JORDAN回车后得到结果如下A110223121用GAUSSJANDAN算得矩阵A的逆矩阵为INVA413513614解法二用通过将矩阵施行初等行变换化成行简化阶梯形的办法,可以借助命令求解,非常方便REFC输入矩阵A110223121A110223121CA,EYELENGTHAC110100223010121001INVARREFCINVA100413010513001614后三列即为其逆阵,检验其正确性CINVA,46C413513614DCAD100010001可见求解正确例24用分解LU的方法求解方程组12342469958107X解解线性方程组中分解的L,U可以实现矩阵A的三角分解,使得ALULUL,U应该是下三角和上三角矩阵的,这样才利于回代求根但是MATLAB中的LU分解与解线性方程组中的L,U不一样MATLAB的LU分解命令调用格式为L,ULUAMATLAB计算出来的L是“准下三角“交换L的行后才能成为真正的下三角阵,U为上三角矩阵,但它们还是满足ALU的先录入矩阵系数A242649615269186151840A242649615269186151840B9232247B9232247将A作分解,方法是使用矩阵分解的LU命令即可LUL,DLUAL03333100000666710000066671000000033331000010000010000000U600001500001800004000000100006000011666700300007000000003333再检验其正确性CLUC242649615269186151840解方程组LYBYLBY470000833332000003333解方程组得到方程组的最终解UXYXUYX05000200003000010000故方程组的最终解为05,231TX例25解方程组123661440X试用平方根法,改进的平方根法和追赶法分别解之解先输入矩阵A6101410114A6101410114B624322B624322平方根法先对A矩阵作分解CHOLESKYLCHOLAL24495040820019579051080037066检验其正确性LLANS60000100000100004000010000010000140000将L转化为下三角矩阵LLL24495000408219579000510837066解方程组LYBYLBY24495117473852526再解方程组,得到最终解TLXYXLYX1000000000230000故平方根法解上述方程的结果为1,023TX改进的平方根法先对矩阵A作LDL分解L,DLDLAL10000000166710000000260910000D6000000038333000137391检验其分解正确性LDLANS6101410114解方程组LXYYLBY623316解方程组TDLXYXDLYX1023故改进的平方根法解上述方程的结果亦为1,023TX追赶法编制追赶法求解该方程的程序如下PURSUEM三对角线性方程组的追赶法解方程组例25输入矩阵CLCA6101410114F624322N,MSIZEA分别取对角元素AZEROS1,NA2NDIAGA,1CDIAGA,1此处用变量D存储A主对角线上的元素,因已用变量B存储方程右边的系数BDIAGAIFB10ERROR主对角元素不能为0RETURNEND初始计算,式231ALPHA1B1BETA1C1/B1按照公式231计算FORI2N1ALPHAIBIAIBETAI1IFALPHAI0ERROR错误在解方程过程中为0RETURNENDBETAICI/ALPHAIEND对最后一行作计算ALPHANBNANBETAN1IFALPHAN0ERROR错误在解方程过程中最后一个为0RETURNEND以下按照公式232计算,解LYFY1F1/B1FORI2NYIFIAIYI1/ALPHAIEND以下按照公式233计算,解UXYXNYNFORIN111XIYIBETAIXI1ENDDISPXDISPX在MATLAB命令窗口输入PURSUE计算结果如下PURSUEA6101410114F624322X1023其中A为系数矩阵,F为矩阵右端的系数,最后计算结果为X由以上计算可知追赶法解该方程的结果亦为1,023TX例26,求1,053TX12,X解输入矩阵X10503X100000500003000计算其1范数NORM_1NORMX,1NORM_118000计算其2范数NORM_2NORMXNORM_211576计算其无穷大范数NORM_INFNORMX,INFNORM_INF1还可以计算其无穷小范数即各分量绝对值中的最小值NORM_MINUSINFNORMX,INFNORM_MINUSINF03000例27210求12,AXA解先输入矩阵A210121012A210121012求A的1范数列和范数NORM_1NORMA,1NORM_14求解A的无穷大范数行和范数NORM_INFNORMA,INFNORM_INF4求A的2范数最大特征值TANORM_2NORMA,2NORM_234142还可以求解A的F范数NORM_FNORMA,FRONORM_F4谱半径可以通过按其概念进行计算对其特征值的绝对值取最大值即可R_AMAXABSEIGAR_A34142例28N阶HILBERT矩阵112341135211NNHNN考查其条件数解上述矩阵为抽象矩阵,而计算机只能进行有限次迭代,我们从考查3,10N其条件数,为此编制如下程序HILB_COND10MHILB_COND10M考查从3阶至10的矩阵2条件数HILBERTFORI310HHILBICONDAICONDH,2ENDDISPNCONDFORI310SSPRINTFDF,I,CONDAIDISPSEND运行后得到如下结果NCOND35240567784155137387395476607250242614951058641005747536735637039381525757527077236494931539864662709401016025391750078617000结果中左边为的阶数,右边为对应的条件数,从以上结果HILBERT2NCODH也可分析可知当N的阶数稍高时,矩阵即出现严重的病态ILBERT例29求213H的条件数2122,CONDHCOND解建立2阶矩阵ILBERTHHILB2H10000050000500003333求其2条件数COND_2CONDH,2COND_2192815求其1条件数COND_1CONDH,1COND_1270000求其无穷大条件数COND_INFCONDH,INFCOND_INF270000还可求其F条件数COND_FCONDH,FROCOND_F19333322MATLAB解线性方程组常用命令介绍1求秩命令RANK在解线性方程组时,为了判断是否存在解,我们应先判断矩阵的秩调用格式为CRANKA2矩阵零空间向量NULL当线性方程组的秩时,方程组有无穷多个解,使用XNULLA可0AXRN得到满足的一个解向量,可为符号矩阵或数值矩阵XNULLA或XNULLSYMA3方阵的LU分解命令LU调用格式为L,ULUAL为准下三角矩阵,U为上三角矩阵,满足ALUL,U,PLUAL为下三角矩阵,U为上三角矩阵,P为变换方阵元素位置的换位阵,满足PAPL其它调用格式请用HELPLU获得更多信息4CHOLSKY分解CHOLLCHOLA其中L为一个下三角矩阵,满足必须为正定矩阵TAL5条件数CONDCCONDA,PA可为向量或矩阵,P取值为下列之一1向量或矩阵的返回1条件数2返回向量或矩阵的2条件数INF返回向量或矩阵的条件数FRO返回向量或矩阵的F条件数6奇异值分解SVDU,S,VSVDA将A分解为正交矩阵U,对角矩阵S和正交矩阵V的乘积,使得AUSVT7线性方程组的符号解LINSOLVEXLINSOLVEA,B它等价于XSYMASYMB返回方程组的符号解3线性方程组的迭代解法31例题解答例31用JACOBI迭代法解以下方程1230150X解编制迭代计算的M文件程序如下JACOBI迭代法求解例31A为方程组的增广矩阵CLCA1021321011512510MAXTIME50最多进行50次迭代EPS1E5迭代误差N,MSIZEAXZEROSN,1迭代初值YZEROSN,1K0进入迭代计算DISP迭代过程X的值情况如下DISPXWHILE1DISPXFORI11NS00FORJ11NIFJISSAI,JXJENDYIAI,N1S/AI,IENDENDFORI11NMAXEPSMAX0,ABSXIYI检查是否满足迭代精度要求ENDIFMAXEPSMAXTIME超过最大迭代次数退出ERROR超过最大迭代次数,退出RETURNENDEND运行该程序结果如下A1021321011512510迭代过程X的值情况如下X000030001500020000080001760026600091801926028640097161970029540098941989729823099621996129938099861998629977099951999529992099981999829997099991999929999100002000030000100002000030000容易看出迭代计算最后结果为1,23TX例32用GAUSSSEIDEL迭代法计算例31并作比较解编制求解程序GAUSS_SEIDELM如下GAUSS_SEIDELMGAUSS_SEIDEL迭代法求解例32A为方程组的增广矩阵CLCFORMATLONGA1021321011512510N,MSIZEA最多进行50次迭代MAXTIME50控制误差EPS10E5初始迭代值XZEROS1,NDISPX迭代次数小于最大迭代次数,进入迭代FORK1MAXTIMEDISPXFORI1NS00FORJ1NIFIJSSAI,JXJ计算和ENDENDXIAI,N1S/AI,I求出此时迭代的值END因为方程的精确解为整数,所以这里将迭代结果向整数靠近的误差作为判断迭代是否停止的条件IFSUMXFLOORX2GAUSS_SEIDELA1021321011512510X000030000000000000015600000000000002684000000000000088040000000000019444800000000002953872000000000098428320000000019922438400000002993754176000000099782418560000019989402547200002999140939008000099970214484480019998545228697602999882238116864099995912838563819999800494888142999983845472654099999439444502819999972634362712999997784263514099999923111360619999996246490722999999696082350099999989453804919999999485158452999999958313948099999998553456419999999929383082999999994282236099999999801588519999999990314012999999999215737099999999972785419999999998671452999999999892429099999999996267219999999999817782999999999985246099999999999488019999999999975012999999999997976099999999999929819999999999996572999999999999722099999999999990419999999999999532999999999999962099999999999998719999999999999942999999999999995099999999999999819999999999999992999999999999999100000000000000020000000000000003000000000000000迭代结果X123可见对此题GAUSS_SEIDEL法的收敛速度还是很快的例33取,用GAUSSSEIDEL迭代法和SOR法计算下列方程组的解4512340XX解GAUSSSEIDEL迭代法可利用上题中的程序,把输入矩阵A换掉就可以了,以下编制求解该程序的SORMSOR法求解例33W145方程组系数矩阵CLCA421242123方程组右端系数B0,2,3W145最大迭代次数MAXTIME100精度要求EPS1E5以15位小数显示FORMATLONGNLENGTHAK0初始迭代值XONESN,1YXDISP迭代过程DISPXWHILE1YXDISPX计算过程FORI1NSBIFORJ1NIFJISSAI,JXJENDENDIFABSAI,IMAXTIMEERROR已达最大迭代次数或矩阵系数近似为0,无法进行迭代RETURNENDSS/AI,IXI1WXIWSENDIFNORMYX,INFGAUSS_SEIDELA421024221233X111075000000000000003750000000000001500000000000000056250000000000005312500000000001541666666666667065104166666666705963541666666671614583333333333070182291666666706582031250000001672743055555556074728732638888907100151909722221722439236111111078561740451388907540283203125001764558015046297081815366391782407913558394820601800288447627315084575003164785908230192396375871830596170307677086915866239571308498774163516951856304498366368088901483276743908726596655669031878111387967082090585767977522208919845338711521896608915839176092014449589537009083767058672731912299302543305093226317856946309222812405563841925608553227410094254275858504409340756559062271936898023465833095126233381957209440801786427021946474230368326095865864691343309525664386408791954597174731730096493251300337209597648438675511961487400246158097025427199531509658708361207371967331981412263097476841341343409710501974128481972289602746377097859749939301809754435510696981976494867177471098184549232921709791701797533441980061950611968098460057752966409823312640708161983087701890432098693755750801609850126296992241985654272302155098891988292515109872870776136531987831346050819099060137531953109892163606851751989678032229960099202768840007809908528603150191991244469676705099323754757668609922410086266961992573188276692099426380138252109934184948296071993700263680578099513431333494809944172885077631994656296783491099587271844975409952645076166231995467244561000099649906494856109959831547547811996155124819374099703035858223409965927417008041996738613994614099748102434905609971098191718351997