1、数值分析实验教程何光辉 重庆大学数学与统计学院二 O 一 O 年三月2目 录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 非
2、线性方程求根的命令 796 插值法 796.1 例题解答 .796.2 Matlab 插值函数介绍 .957 数据拟合和最佳平方逼近 .967.1 例题解答 .967.2 Matlab 数据拟合命令介绍 1068 数值积分与数值微分 .1078.1 例题解答 1079 常微分方程数值解法 .1309.1 例题解答 .1309.2 Matlab 常微分方程数值解常用命令介绍 14630 Matlab 介绍入门知识1. Matlab简介MATLAB的含义是矩阵实验室(MATRIX LABORATORY),主要用于方便矩阵的存取,其基本元素是无须定义维数的矩阵.MATLAB自问世以来,就是以数值计算
3、称.MATLAB进行数值计算的基本单位是复数数组(或称阵列),这使得MATLAB高度“向量化”.经过十几年的完善和扩充,现已发展成为线性代数课程的标准工具.由于它不需定义数组的维数,并给出矩阵函数、特殊矩阵专门的库函数,使之在求解诸如信号处理、建模、系统识别、控制、优化等领域的问题时,显得大为简捷、高效、方便,这是其它高级语言所不能比拟的.MATLAB中包括了被称作工具箱(TOOLBOX)的各类应用问题的求解工具.工具箱实际上是对MATLAB进行扩展应用的一系列MATLAB函数(称为M文件),它可用来求解各类学科的问题,包括信号处理、图象处理、控制系统辨识、神经网络等.随着MATLAB版本的不
4、断升级,其所含的工具箱的功能也越来越丰富,因此,应用范围也越来越广泛,MATLAB提供的工具箱已覆盖信号处理、系统控制、统计计算、优化计算、神经网络、小波分析、偏微分方程、模糊逻辑、动态系统模拟、系统辨识和符号运算等领域.当前它的使用范围涵盖了工业、电子、医疗、建筑等各行各业.MATLAB中包括了图形界面编辑GUI,让使用者也可以象VB、VC、VJ、DELPHI等那样进行一般的可视化的程序编辑.在命令窗口(matlab command window)键入simulink,就出现(SIMULINK) 窗口.以往十分困难的系统仿真问题,用SIMULINK只需拖动鼠标即可轻而易举地解决问题,这也是近
5、来受到重视的原因所在.MATLAB 语言由美国 The MathWorks 开发,最早是由C.Moler用Fortran语言编写的,用来方便地调用LINPACK和EISPACK矩阵代数软件包的程序.后来他创立了MATHHWORKS公司,对MATLAB作了大量的、坚持不懈的改进.Cleve B.Moler是The MathWork公司的主席和首席科学家.曾任密歇系教授.他在两个计算机硬件制造商Intel公司的Hypercube组织和Arden Computers 公司工作了五年.他的主要专业兴趣在于数值分析和科学计算.他是MATLAB软件的创始者,也是著名的矩阵计算软件包LINPACK和EISP
6、ACK的著作这一,已撰写了三本有相关数值方法的教材.同时,他在SIAM(美国工业与应用数学学会)历任期刊编辑、委员会成员和副总裁,并从1996年开始担任理事会成员.2. Matlab 入门知识Matlab 变量名是以字母开头,后接字母、数字或下划线的字符序列,最多 63 个字符.在MATLAB 中,变量名区分字母的大小写.赋值语句:变量=表达式 或 表达式其中表达式是用运算符将有关运算量连接起来的式子,其结果是一个矩阵.clear命令用于删除MATLAB工作空间中的变量.who和whos这两个命令用于显示在MATLAB工作空间中已经驻留的变量名清单.who命令只显示出驻留变量的名称,whos在
7、给出变量名的同时,还给出它们的大小、所占字节数及数据类型等信息.利用MAT文件可以把当前MATLAB工作空间中的一些有用变量长久地保留下来,扩展名是.mat.MAT文件的生成和装入由save和load命令来完成.常用格式为:save 文件名 变量名表 -append-asciiload 文件名 变量名表 -ascii其中,文件名可以带路径,但不需带扩展名.mat,命令隐含一定对.mat文件进行操作.变量名表中的变量个数不限,只要内存或文件中存在即可,变量名之间以空格分隔.当变量名表省4略时,保存或装入全部变量.-ascii选项使文件以ASCII格式处理,省略该选项时文件将以二进制格式处理.sa
8、ve命令中的-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 =1234若不想输出结果,在每一条语句后用分号作为结束符,若留空或用逗号结束,则在执行该语句后会把结果输出来. a+b; a+bans =
9、2 3 4 5 6(2) 矩阵的创建直接输入:最简单的建立矩阵的方法是从键盘直接输入矩阵的元素.具体方法如下:将矩阵的元素用方括号括起来,按矩阵行的顺序输入各元素,同一行的各元素之间用空格或逗号分隔,不同行的元素之间用分号分隔. A=1 2 3;4 5 6;2 3 5A =1 2 34 5 62 3 5利用矩阵函数创建: B=magic(3)%魔方阵B =8 1 63 5 74 9 2 C=hilb(3)%3 阶 Hilbert 矩阵C =51.0000 0.5000 0.33330.5000 0.3333 0.25000.3333 0.2500 0.2000Matlab 中用%引导注释其它创
10、建矩阵函数还有: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): 生
11、成一个以x为第一列,y为第一行的托普利兹矩阵compan(p): 生成伴随矩阵, p是一个多项式的系数向量,高次幂系数排在前,低次幂排在后.pascal(n): 生成一个n阶帕斯卡矩阵.compan: 生成伴随矩阵(3) 矩阵运算MATLAB 的基本算术运算有:(加)、(减)、*(乘)、/(右除)、(左除)、(乘方).加法: A+Bans =9 3 97 10 136 12 7减法: B-Aans =7 -1 3-1 0 12 6 -3乘法: A*Bans =26 38 2671 83 7145 62 43除法: magic(3)/hilb(3)ans =61.0e+003 *0.2160 -
12、1.1760 1.14000.0570 -0.4080 0.4500-0.2280 1.2240 -1.1400在MATLAB中,有一种特殊的运算,因为其运算符是在有关算术运算符前面加点,所以叫点运算.点运算符有.*、./、.和.两矩阵进行点运算是指它们的对应元素进行相关运算,要求两矩阵的维参数相同. A.*Bans =8 2 1812 25 428 27 10MATLAB提供了6种关系运算符:(大于)、=(大于或等于)、=(等于)、=(不等于). ABans =0 1 01 0 00 0 1MATLAB提供了3种逻辑运算符: plot(x,peaks); box on; title(绘制混合
13、图形); xlabel(X 轴); ylabel(Y 轴);绘制图像为:0 5 10 15 20 25 30 35 40 45 50-8-6-4-20246810 一一一一一一X一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);10-6 -4 -2 0 2 4 6-6-4-202
14、46(5) 二维符号函数曲线专用命令 ezplotf = f(x)时 :ezplot(f):在默认区间-2 figure;ezplot(cos(tan(pi*x), 0,1);110 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1-0.500.51xcos(tan( x)(6) 图形窗口的分割 subplotsubplot(m,n,p)该函数将当前图形窗口分成 mn 个绘图区,即每行 n 个,共 m 行,区号按行优先编号,且选定第 p 个区为当前活动区.在每一个绘图区允许以不同的坐标系单独绘制图形.(7) 其他坐标系下的二维数据曲线图对数坐标图形:semilog
15、x(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,z
16、n,选项 n) 其中每一组 x,y,z 组成一组曲线的坐标参数,选项的定义和 plot 函数相同.当 x,y,z 是同维向量时,则 x,y,z 对应元素构成一条三维曲线.当 x,y,z 是同维矩阵时,则以 x,y,z 对应列元素绘制三维曲线,曲线条数等于矩阵列数.12 t=0:0.1:8*pi; plot3(sin(t),cos(t),t);-1-0.500.51-1-0.500.51051015202530(9) 产生三维数据在 MATLAB 中,利用 meshgrid 函数产生平面区域内的网格坐标矩阵.其格式为:X,Y=meshgrid(x,y);语句执行后,矩阵 X 的每一行都是向量 x
17、,行数等于向量 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
18、) 其他三维绘图指令介绍bar3 函数绘制三维条形图,常用格式为bar3(y)bar3(x,y)13stem3 函数绘制离散序列数据的三维杆图,常用格式为: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
19、)程序的暂停: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
20、语句组 2elseif 条件 m语句组 melse语句组 nend语句用于实现多分支选择结构.14(7)switch 语句:switch 表达式case 表达式 1语句组 1case 表达式 2语句组 2case 表达式 m语句组 motherwise语句组 nend(8)try 语句语句格式为:try语句组 1catch语句组 2endtry 语句先试探性执行语句组 1,如果语句组 1 在执行过程中出现错误,则将错误信息赋给保留的 lasterr 变量,并转去执行语句组 2.(9)for 语句for 语句的格式为:for 循环变量=表达式 1:表达式 2:表达式 3循环体语句end其中表达式
21、 1 的值为循环变量的初值,表达式 2 的值为步长,表达式 3 的值为循环变量的终值.步长为 1 时,表达式 2 可以省略.for 语句更一般的格式为:for 循环变量=矩阵表达式循环体语句end执行过程是依次将矩阵的各列元素赋给循环变量,然后执行循环体语句,直至各列元素处理完毕.(10)while 语句while 语句的一般格式为:while (条件)循环体语句end其执行过程为:若条件成立,则执行循环体语句,执行后再判断条件是否成立,如果不成立则跳出循环.(11)break 语句和 continue 语句与循环结构相关的语句还有 break 语句和 continue 语句.它们一般与 if
22、 语句配合使用.break 语句用于终止循环的执行.当在循环体内执行到该语句时,程序将跳出循环,继15续执行循环语句的下一语句.continue 语句控制跳过循环体中的某些语句 .当在循环体内执行到该语句时,程序将跳过循环体中所有剩下的语句,继续下一次循环.(12)循环的嵌套如果一个循环结构的循环体又包括一个循环结构,就称为循环的嵌套,或称为多重循环结构.(13)函数文件的基本结构函数文件由 function 语句引导,其基本结构为function 输出形参表 =函数名(输入形参表)注释说明部分函数体语句其中以 function 开头的一行为引导行,表示该 M 文件是一个函数文件.函数名的命名
23、规则与变量名相同.输入形参为函数的输入参数,输出形参为函数的输出参数.当输出形参多于一个时,则应该用方括号括起来.(14)函数调用函数调用的一般格式是:输出实参表=函数名(输入实参表)注意的是,函数调用时各实参出现的顺序、个数,应与函数定义时形参的顺序、个数一致,否则会出错.函数调用时,先将实参传递给相应的形参,从而实现参数传递,然后再执行函数的功能.在 MATLAB 中,函数可以嵌套调用,即一个函数可以调用别的函数,甚至调用它自身.一个函数调用它自身称为函数的递归调用.(15)函数参数的可调性在调用函数时,MATLAB 用两个永久变量 nargin 和 nargout 分别记录调用该函数时的
24、输入实参和输出实参的个数.只要在函数文件中包含这两个变量,就可以准确地知道该函数文件被调用时的输入输出参数个数,从而决定函数如何进行处理.(16)全局变量与局部变量全局变量用 global 命令定义,格式为:global 变量名(17)程序调试Debug 菜单项:该菜单项用于程序调试,需要与 Breakpoints 菜单项配合使用.Breakpoints 菜单项:该菜单项共有 6 个菜单命令,前两个是用于在程序中设置和清除断点的,后 4 个是设置停止条件的,用于临时停止 M 文件的执行,并给用户一个检查局部变量的机会,相当于在 M 文件指定的行号前加入了一个 keyboard 命令.调试命令:
25、除了采用调试器调试程序外,MATLAB 还提供了一些命令用于程序调试.命令的功能和调试器菜单命令类似,具体使用方法请读者查询 MATLAB 帮助文档.161 绪论1.1 例题解答例 1.1 计算 , .sinx04解:创建符号函数: 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求泰勒级数在 处的函数值:0.5x subs(h,x,0.5) ans =0.479425533234127也可以通过内联函数来求解:H=inline(h) H =Inline func
26、tion:H(x) = x-1./6.*x.3+1./120.*x.5-1./5040.*x.7 feval(H,0.5) ans =0.479425533234127例 1.2 计算积分值 .10Idx解:解法一:( 符号法) : 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:5I=I+(f(i)+f(i+1)/2*0.2; %将每一小曲边的梯形累加起来作为积分值End17 vpa(I,9) %取结果的小数精度为 9 位小数ans
27、=.695634921例 1.3 略例 1.4 不用开平方根计算 的值.(0)a解:解法一(符号法): A=sym(a); sqrt(A)ans =a(1/2)解法二(数值法):按以下迭代公式迭代计算近似值: 1(),01,2.2kkaxx建立函数文件 msqrt.mfunction x=msqrt(x0,a)%用迭代法近似计算平方根%x0 为初始迭代值,a 为开平方数format long;x=zeros(20,1);x(1)=x0;for i=2:20x(i)=1/2*(x(i-1)+a/x(i-1);enddisp(x);用编写的函数计算 , :302x msqrt(2,3);2.000
28、0000000000001.7500000000000001.7321428571428571.7320508100147271.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.732050807568877181.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.732050807568
29、8771.7320508075688771.732050807568877上述结果为迭代过程计算的中间结果,分析数据可知迭代收敛速度快,只需四次计算即可计算出较为准确的数值.例 1.5 略.例 1.6 计算 ,视已知数为精确数,用 4 位浮点数计算 .175960解:直接在 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=sy
30、m(1/759),b=sym(1/760),c=a-ba =1/759b =1/760c =1/576840可见符号运算准确但耗费运算时间.例 1.7 略.例 1.8 解方程 .2180x解:符号法解方程: x=solve(x2-18*x+1,x)x =199+4*5(1/2)9-4*5(1/2)将结果保留小数点 6 位: vpa(x,6)ans =17.9443.5572e-11.2 Matlab 中数值计算精度1. Matlab 中有三种运算精度 ,它们分别为数值算法、符号算法和可控精度算法,将它们分别介绍如下:(1) 数值算法把每个数取为 16 位, 计算按浮点运算进行,它是运算速度最快
31、的一种算法.(2) 符号算法把每个数都变为符号量,运算按有理量计算进行,它的优点是能够得到精确结果, 缺点是占用空间大, 并且运算速度最慢 .(3) 可控精度算法介于上述两种算法之间,它能够使运算在可控的精度下进行计算.2. Matlab 的数据显示格式 ,列表如下 :表 Matlab 数据显示格式命令命令 意义 举例( )format short 短格式方式,显示 5 位定点十进制数 3.1416format long 长格式方式,显示 15 位定点十进制数 3.141592653589793format short e 最优化短格式显示,5 位加指数 3.1416e+000format l
32、ong e 最优格式,15 位加指数 3.141592653589793e+000format short g 5 位定点或浮点格式 3.1416format long g 对双精度,显示 15 位定点或浮点格式,对单精度,显示 7 位定点或浮点格式3.14159265358979format short eng 至少 5 位加 3 位指数 3.1416e+000format long eng 16 位加至少 3 位指数 3.14159265358979e+000format hex 十六进制格式方式 400921fb54442d18format bank 银行格式.按元、角、分(小数点后具有
33、两位)的固定格式3.14format + +格式,以, 和空格分别表示中的正数,负数和零元素+format 缺省时为默认短格式方式与 format short 相同3.1416format rat 分数格式形式.用有理数逼近显示数据 355/113 format loose 松散格式.数据之间有空行format compact 紧凑格式.数据之间无空行vpa(date,n) 将数据 date 以 n 位有效数字显示 vpa(pi,5)= 3.141620format 并不影响 matlab 如何计算和存储变量的值.对浮点型变量的计算,即单精度或双精度,按合适的浮点精度进行,而不论变量是如何显示
34、的.对整型变量采用整型数据.整型变量总是根据不同的类(class)以合适的数据位显示.3. Matlab 的特殊变量ans :对最近输入的反应 computer :当前计算机类型 eps :浮点精度 flops :计算浮点操作次数,现已不再常用 i :虚部单位 inf :无穷大 inputname :输入参数名 j :虚部单位 nan :非数值 nargin :输入参数的数目 nargout :输出参数的数目(用户定义函数) pi :圆周率 realmax :最大正浮点数 realmin :最小正浮点数 varargin varargout :返回参数数目(matlab函数) cputime:
35、CPU工作时间212 线性方程组的直接解法2.1 例题解答例 2.1 用 Gauss 消元法解方程组 : 123175649xx解:直接建立求解该方程组的 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=nerror(矩阵 A 的行数和列数必须相同);return;endif m=size(b)error(b 的大小必须和 A 的行数或 A 的列数相同);return
36、;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%回代结果22x=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%显示计
37、算结果disp(x=);disp(x);直接运行上面的 M 文件或在 Matlab 命令窗口中直接输入 Gauss 即可得出结果.在 Matlab 命令窗口中输入 Gauss 得出结果如下: Gaussx=2.00001.0000-1.0000扩展:Matlab 求解线性方程的几种命令如下(方程组的一般形式可用矩阵和向量表示成: ,但运用下列方法的前提必须保证所求解的方程为恰定方程,即方程组存在唯一的一Axb组解):运用求逆思想 : 或 ;1Ab()*xinvb(1)*xAb左除法,原理上是运用高斯消元法求解, 但 Matlab 在实际执行过程中是通过 分解法进行LU的(即先将矩阵 A 作 分
38、解,再回代计算): ;LU用符号矩阵法计算,这种计算方法最接近精确值, 但计算速度最慢: ()(xsymAb通过将矩阵施行初等行变换化成行简化阶梯形的办法,可以这样实现之: ; ,C.()refC上面四种常用的办法示例如下: A=1 2 3;2 7 5;1 4 9 %上面示例方程组系数A =1 2 32 7 51 4 9 b=1 6 -3 %方程组右端的系数b =16-3 x1_1=inv(A)*b,x1_2=A(-1)*b %方法一,求逆思想x1_1 =232.00001.0000-1.0000x1_2 =2.00001.0000-1.0000 x2=Ab %方法二,左除思想x2 =21-1
39、 x3=sym(A)sym(b) %方法三, 符号法x3 =21-1 C=A,b,rref(C) %方法四, 行简化阶梯形思想, 最后输出结果的一列为解C =1 2 3 12 7 5 61 4 9 -3ans =1 0 0 20 1 0 10 0 1 -1例 2.2 用选列主元的 Gauss 消元法解如下方程 :120013x解:直接建立求解该方程的 M 文件 Gauss_line.m,求解程序编制如下 :%求解例题 2.2%高斯列主元消元法求解线性方程组 Ax=b%A 为输入矩阵系数,b 为方程组右端系数%方程组的解保存在 x 变量中format long;%设置为长格式显示,显示 15 位
40、小数A=0.00001,1;2,1b=1.00001,324m,n=size(A);%先检查系数正确性if m=nerror(矩阵 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
41、,k)=0) if(m=k)A(k m,:)=A(m k,:); %换行endA(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:1x(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.000010000000000 1.000000000
42、0000002.000000000000000 1.000000000000000b =1.0000100000000003.000000000000000X=251.0000000000000001.000000000000000例 2.3 用 消元法思想求GausJordn1023A的逆阵.解:解法一:直接建立求解的 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=merror(A 必须为方阵);return;endA(:,n
43、+1:2*n)=eye(n);%构造增广矩阵for k=1:nl,m=max(abs(A(k:n,k);%按列选主元if A(k+m-1,k)=0error(找到列最大的元素为零, 错误);return;endif m=1 %交换Temp=A(k,:);A(k,:)=A(k+m-1,:);A(k+m-1,:)=Temp;endfor i=1:nif i=kA(i,:)=A(i,:)-A(k,:)*A(i,k)/A(k,k);endend endfor i=n:(-1):1A(i,:)=A(i,:)/A(i,i);endA(:,1:n)=;disp(A=);26disp(A1);disp(用 G
44、auss-Jandan 算得矩阵 A 的逆矩阵为:);disp(inv(A)=);disp(A); clear Temp i k l m n;%清除临时变量在 Matlab 命令窗口中输入 Gauss_Jordan 回车后得到结果如下:A=1 -1 02 2 3-1 2 1用 Gauss-Jandan 算得矩阵 A 的逆矩阵为:inv(A)=-4 1 -3-5 1 -36 -1 4解法二:用通过将矩阵施行初等行变换化成行简化阶梯形的办法,可以借助 命令求解,非()refC常方便:输入矩阵: A=1 -1 0;2 2 3;-1 2 1A =1 -1 02 2 3-1 2 1 C=A,eye(le
45、ngth(A)C =1 -1 0 1 0 02 2 3 0 1 0-1 2 1 0 0 1 invA=rref(C)invA =1 0 0 -4 1 -30 1 0 -5 1 -30 0 1 6 -1 4后三列即为其逆阵,检验其正确性: c=invA(:,4:6)c =-4 1 -3-5 1 -36 -1 4 d=c*Ad =271 0 00 1 00 0 1可见求解正确.例 2.4 用分解 LU 的方法求解方程组 12342469958107x解: 解线性方程组中 分解的 L,U 可以实现矩阵 A 的三角分解,使得 A=L*U. L,U 应该是下三LU角和上三角矩阵的,这样才利于回代求根.但
46、是 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 64 9 6 152 6 9 186 15 18 40 b=9 23 22 47b =9232247将 A 作 分解,方法是使用矩阵分解的 LU 命令即可:LUL,D=lu(A)L =0.3333 1.0000 -0.
47、6667 1.00000.6667 1.0000 0 00.3333 -1.0000 1.0000 01.0000 0 0 0U =6.0000 15.0000 18.0000 40.0000280 -1.0000 -6.0000 -11.66670 0 -3.0000 -7.00000 0 0 -0.3333再检验其正确性:C=L*U C =2 4 2 64 9 6 152 6 9 186 15 18 40解方程组 Lyby=Lb y =47.0000-8.3333-2.00000.3333 解方程组 得到方程组的最终解:Uxyx=Uy x =0.50002.00003.0000-1.0000 故方程组的最终解为: (0.5,231)Tx例 2.5 解方程组 123606144x试用平方根法,改进的平方根法和追赶法分别解之.解:先输入矩阵:A=6 1 0;1 4 1;0 1 14 A =6 1 01 4 10 1 14 b=6 24 322 b =29624322平方根法:先对 A 矩阵作 分解:CholeskyL=chol(A) L =2.4495 0.4082 00 1.9579 0.51080 0 3.7066 检验其正确性L*L ans =6.