数值分析教案 - 图文

数值分析实验

董海云

数理学院数学实验教学中心

2008年3月

目 录

0 Matlab介绍入门知识 ........................................................................................................... 3 1 绪论 ..................................................................................................................................... 17

1.1 例题解答 ............................................................................................................... 17 1.2 Matlab中数值计算精度 ...................................................................................... 20 2 线性方程组的直接解法 ..................................................................................................... 22

2.1 例题解答 ............................................................................................................... 22 2.2 Matlab解线性方程组常用命令介绍 .................................................................. 36 3 线性方程组的迭代解法 ..................................................................................................... 37

3.1 例题解答 ................................................................................................................ 37 3.2 Matlab迭代解法用到的函数介绍 ...................................................................... 53 4 方阵特征值和特征向量的计算 ......................................................................................... 54

4.1 例题解答 ................................................................................................................ 54 4.2 Matlab关于方阵特征值为特征向量函数介绍 .................................................. 62 5 非线性方程求根 ................................................................................................................. 63

5.1 例题解答 ................................................................................................................ 63 5.2 Matlab非线性方程求根的命令 .......................................................................... 84 6 插值法 ............................................................................................................................... 85

6.1 例题解答 ................................................................................................................ 85 6.2 Matlab插值函数介绍 ......................................................................................... 100 7 数据拟合和最佳平方逼近 ............................................................................................. 102

7.1 例题解答 .............................................................................................................. 102 7.2 Matlab数据拟合命令介绍 ................................................................................ 112 8 数值积分与数值微分 ..................................................................................................... 113

8.1 例题解答 ............................................................................................................. 113 9 常微分方程数值解法 ..................................................................................................... 136

9.1 例题解答 .............................................................................................................. 136 9.2 Matlab常微分方程数值解常用命令介绍 ........................................................ 153

2

0 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工作空间中的一些有用变量长久地保留下

3

来，扩展名是.mat.MAT文件的生成和装入由save和load命令来完成.常用格式为：

save 文件名 [变量名表] [-append][-ascii] load 文件名 [变量名表] [-ascii]

其中，文件名可以带路径，但不需带扩展名.mat，命令隐含一定对.mat文件进行操作.变量名表中的变量个数不限，只要内存或文件中存在即可，变量名之间以空格分隔.当变量名表省略时，保存或装入全部变量.-ascii选项使文件以ASCII格式处理，省略该选项时文件将以二进制格式处理.save命令中的-append选项控制将变量追加到MAT文件中. (1) 向量的创建 用步长生成法：

数组=初值:步长(增量):终值 >> a=1:0.5:3 a =

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+b ans =

2 3 4 5 6 (2) 矩阵的创建 直接输入:

最简单的建立矩阵的方法是从键盘直接输入矩阵的元素.具体方法如下：将矩阵的元素用方括号括起来，按矩阵行的顺序输入各元素，同一行的各元素之间用空格或逗号分隔，不同行的元素之间用分号分隔. >> A=[1 2 3;4 5 6;2 3 5] A =

1 2 3 4 5 6 2 3 5 利用矩阵函数创建: >> B=magic(3)%魔方阵

4

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.2500

0.3333 0.2500 0.2000 Matlab中用%引导注释 其它创建矩阵函数还有:

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+B ans =

9 3 9 7 10 13 6 12 7 减法： >> B-A ans =

7 -1 3 -1 0 1 2 6 -3 乘法：

5

>> A*B ans =

26 38 26 71 83 71

45 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.*B ans =

8 2 18 12 25 42 8 27 10

MATLAB提供了6种关系运算符：<(小于)、<=(小于或等于)、>(大于)、>=(大于或等于)、==(等于)、～=(不等于).

>> A>B ans =

0 1 0 1 0 0 0 0 1

MATLAB提供了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矩阵第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表示某一维的末尾元素下标.

6

(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*B ans =

1.0000 0 -0.0000 -0.0000 1.0000 0

0.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)

7

计算一般函数值:eval(f) 求虚部函数: imag(x) 求实部函数: real(x) 角相位函数:angle(x) 共轭复数函数:conj(x) 沿零方向取整 :fix (x) 舍入取整:round(x)

沿负无穷大方向取整:floor (x) 沿正无穷大方向取整:ceil(x) 求除法的余数: rem 符号函数:sign(x) 最大公约数:gcd() 4. 图形可视化

(1) 二维绘图指令plot

plot函数的基本调用格式为: 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 < r 虚线 圆圈 朝左三红色 角 X > c ： 点线 叉字符 朝右三青色 角 -. + p m 点划线 加号 五角星 洋红 * h y 星号 六角形 黄色 s k 方块 黑色 d w 菱形 白色

8

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轴'); 绘制图像为:

蓝色 9

绘制混合图形108642Y轴0-2-4-6-80510152025X轴3035404550(4) 二维数值函数的专用绘图函数fplot

fplot(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

6420-2-4-6-6-4-20246(5) 二维符号函数曲线专用命令ezplot

f = f(x)时:

ezplot(f)：在默认区间-2π

ezplot(f)：在默认区间-2π

ezplot(f, [xmin,xmax,ymin,ymax])：在区间xmin

ezplot(f, [a,b])：在区间a

ezplot(x,y)：在默认区间0

ezplot(x,y, [tmin,tmax])：在区间tmin < t < tmax绘制x=x(t)和y=y(t)的图形

>> figure;ezplot('cos(tan(pi*x))',[ 0,1]);

11

cos(tan(? x))10.50-0.5-100.10.20.30.40.5x0.60.70.80.91

(6) 图形窗口的分割subplot

subplot(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) 三维曲线plot3

plot3(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对应列元素绘制三维曲线，曲线条数等于矩阵列数.

12

>> t=0:0.1:8*pi;

>> plot3(sin(t),cos(t),t);

30252015105010.50-0.5-1-1-0.50.501(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函数绘制三维条形图，常用格式为

13

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

14

语句组n end

语句用于实现多分支选择结构. (7)switch语句:

switch 表达式 case 表达式1 语句组1 case 表达式2 语句组2 ……

case 表达式m 语句组m otherwise 语句组n end (8)try语句

语句格式为： try

语句组1 catch

语句组2 end

try语句先试探性执行语句组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

其执行过程为：若条件成立，则执行循环体语句，执行后再判断条件是否成立，如果不成立则跳出循环.

15

(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帮助文档.

16

1 绪论 1.1 例题解答

例1.1 计算sinx,x?[0,?4].

解:

创建符号函数: >> syms x; >>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 求泰勒级数在x?0.5处的函数值： >> subs(h,x,0.5) ans =

0.479425533234127 也可以通过内联函数来求解: >>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.479425533234127

例 1.2 计算积分值I??1101?xdx. 解:

解法一:( 符号法)： >> 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

17

I=I+(f(i)+f(i+1))/2*0.2; %将每一小曲边的梯形累加起来作为积分值 End

>> vpa(I,9) %取结果的小数精度为9位小数 ans =

.695634921 例 1.3略

例 1.4 不用开平方根计算a(a?0)的值. 解:

解法一(符号法): >> A=sym('a'); >> sqrt(A) ans = a^(1/2)

解法二(数值法):

按以下迭代公式迭代计算近似值:

x12(xak?1?k?x),k?0,1,2,...

k建立函数文件msqrt.m

function 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)); end

disp(x);

用编写的函数计算3,x0?2: >> msqrt(2,3);

2.000000000000000 1.750000000000000 1.732142857142857 1.732050810014727 1.732050807568877

18

1.732050807568877 1.732050807568877 1.732050807568877 1.732050807568877 1.732050807568877 1.732050807568877 1.732050807568877 1.732050807568877 1.732050807568877 1.732050807568877 1.732050807568877 1.732050807568877 1.732050807568877 1.732050807568877 1.732050807568877

上述结果为迭代过程计算的中间结果,分析数据可知迭代收敛速度快,只需四次计算即可计算出较为准确的数值.

例 1.5 略.

11?例 1.6 计算,视已知数为精确数,用4位浮点数计算. 759760解:

直接在Matlab中输入式子: >> 1/759-1/760 ans =

1.7336e-006

若先转化为浮点数再运算可得: >> a=1/759,b=1/760,a-b a =

0.0013 b =

0.0013 ans =

1.7336e-006

可见Matlba在计算时,数据结构都取为双精度而提高了运算准确度.若以符号运算计算之,有:

>> a=sym('1/759'),b=sym('1/760'),c=a-b a = 1/759 b = 1/760 c =

1/576840

可见符号运算准确但耗费运算时间.

19

例 1.7 略.

例 1.8 解方程x2?18x?1?0. 解:

符号法解方程:

>> 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-1

1.2 Matlab中数值计算精度

1. Matlab中有三种运算精度,它们分别为数值算法、符号算法和可控精度算法,将它们分别介绍如下:

(1) 数值算法把每个数取为16位,计算按浮点运算进行,它是运算速度最快的一种算法.

(2) 符号算法把每个数都变为符号量,运算按有理量计算进行,它的优点是能够得到精确结果,缺点是占用空间大,并且运算速度最慢.

(3) 可控精度算法介于上述两种算法之间,它能够使运算在可控的精度下进行计算.

2. Matlab的数据显示格式,列表如下:

表 Matlab数据显示格式命令 命令 意义 举例(?) format short 短格式方式，显示5位定点十进制3.1416 数 format long 长格式方式，显示15位定点十进3.141592653589793 制数 format short e 3.1416e+000 最优化短格式显示,5位加指数 format long e 3.141592653589793e+000 最优格式,15位加指数 format short g 3.1416 5位定点或浮点格式 format long g 对双精度，显示15位定点或浮点3.14159265358979 格式，对单精度，显示7位定点或浮点格式 format short eng 至少5位加3位指数 3.1416e+000

20

3.000000000000000 X=

1.000000000000000 1.000000000000000

例 2.3 用Gauss?Jordan消元法思想求

?1?10?? A??223?????121??的逆阵.

解:

解法一:直接建立求解的M文件:Gauss_Jordan.m,源程序如下: %Gauss-Jordan法求例2.3 clc;

A=[1 -1 0;2 2 3;-1 2 1];

A1=A;%先保存原来的方阵A [n,m]=size(A); if n~=m

error('A必须为方阵'); return; end

A(:,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 end

for i=n:(-1):1

A(i,:)=A(i,:)/A(i,i); end

26

A(:,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

解法二:用通过将矩阵施行初等行变换化成行简化阶梯形的办法,可以借助

rref(C)命令求解,非常方便:

输入矩阵:

>> A=[1 -1 0;2 2 3;-1 2 1] A =

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

27

>> d=c*A d =

1 0 0 0 1 0 0 0 1

可见求解正确.

例 2.4 用分解LU的方法求解方程组

?2426??x1??9????x???49615???2???23??26918??x3??22???????6151840???x4??47?

解:

解线性方程组中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 40] A =

2 4 2 6 4 9 6 15 2 6 9 18 6 15 18 40 >> b=[9 23 22 47]' b = 9 23 22 47

将A作LU分解,方法是使用矩阵分解的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 0

28

U =

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 解方程组Ly?b >>y=L\\b y =

47.0000 -8.3333 -2.0000 0.3333

解方程组Ux?y得到方程组的最终解: >>x=U\\y x =

0.5000 2.0000 3.0000 -1.0000

故方程组的最终解为: x?(0.5,2,3,?1)T

例 2.5 解方程组

??610???141??x1???x???6??2???24? ??0114????x3????322??试用平方根法,改进的平方根法和追赶法分别解之.

解:

先输入矩阵:

>>A=[6 1 0;1 4 1;0 1 14] A =

6 1 0 1 4 1 0 1 14

29

>>b=[6 24 322]' b = 6 24 322 平方根法:

先对A矩阵作Cholesky分解:

>>L=chol(A) L =

2.4495 0.4082 0 0 1.9579 0.5108 0 0 3.7066 检验其正确性 >>L'*L ans =

6.0000 1.0000 0 1.0000 4.0000 1.0000 0 1.0000 14.0000 将L转化为下三角矩阵: >>L=L' L =

2.4495 0 0 0.4082 1.9579 0 0 0.5108 3.7066 解方程组Ly?b: >>y=L\\b y =

2.4495 11.7473 85.2526

再解方程组LTx?y,得到最终解: >>x=L'\\y x =

1.0000 -0.0000 23.0000

故平方根法解上述方程的结果为: x?(1,0,23)T改进的平方根法:

先对矩阵A作LDL分解: >> [L,D]=ldl(A)

30

L =

1.0000 0 0 0.1667 1.0000 0 0 0.2609 1.0000 D =

6.0000 0 0 0 3.8333 0 0 0 13.7391 检验其分解正确性: >>L*D*L' ans =

6 1 0 1 4 1 0 1 14 解方程组Lx?y: >> y=L\\b y = 6 23 316

解方程组DLTx?y: >> x=(D*L')\\y x = 1 0 23

故改进的平方根法解上述方程的结果亦为: x?(1,0,23)T

追赶法:

编制追赶法求解该方程的程序如下: %pursue.m

%三对角线性方程组的追赶法解方程组例2.5 %输入矩阵 clc;

A=[6 1 0;1 4 1;0 1 14] f=[6 24 322] [n,m]=size(A);

%分别取对角元素 a=zeros(1,n); a(2:n)=diag(A,-1); c=diag(A,1);

%此处用变量d存储A主对角线上的元素,因已用变量b存储方程右边的系数

31

b=diag(A); if b(1)==0

error('主对角元素不能为0'); return; end

%初始计算,式(2.31) alpha(1)=b(1); beta(1)=c(1)/b(1); %按照公式(2.31)计算 for i=2:n-1

alpha(i)=b(i)-a(i)*beta(i-1); if alpha(i)==0

error('错误:在解方程过程中α为0'); return; end

beta(i)=c(i)/alpha(i); end

%对最后一行作计算

alpha(n)=b(n)-a(n)*beta(n-1); if alpha(n)==0

error('错误:在解方程过程中最后一个α为0'); return; end

%以下按照公式(2.32)计算,解Ly=f y(1)=f(1)/b(1); for i=2:n

y(i)=(f(i)-a(i)*y(i-1))/alpha(i); end

%以下按照公式(2.33)计算,解Ux=y X(n)=y(n); for i=n-1:-1:1

X(i)=y(i)-beta(i)*X(i+1); end

disp('X='); disp(X);

在Matlab命令窗口输入pursue计算结果如下: >>pursue A =

6 1 0 1 4 1 0 1 14 f =

6 24 322

32

X=

1 0 23

其中A为系数矩阵,f为矩阵右端的系数,最后计算结果为X. 由以上计算可知追赶法解该方程的结果亦为: x?(1,0,23)T

例 2.6 x?(1,0.5,?0.3)T,求x1,x2,x?.

解:

输入矩阵:

>>x=[1 0.5 -0.3] x =

1.0000 0.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

???210??1?21?? ??01?2??求A1,x?,A2,?(A).

解:

先输入矩阵:

>>A=[-2 1 0;1 -2 1;0 1 -2] A =

-2 1 0 1 -2 1 0 1 -2

33

求A的1-范数(列和范数): >>norm_1=norm(A,1) norm_1 = 4

求解A的无穷大范数(行和范数): >>norm_inf=norm(A,inf) norm_inf = 4

求A的2-范数(ATA最大特征值): >>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.8 n阶Hilbert矩阵

??1??1?2?Hn??1?3????1??n考查其条件数. 解:

121314?1n?1131415?1n?2?????1?n??1?n?1??1? n?2????1??2n?1?上述矩阵为抽象矩阵,而计算机只能进行有限次迭代,我们从n?3,?,10考查其条件数,为此编制如下程序Hilb_cond10.m:

34

% Hilb_cond10.m

%考查从3阶至10的Hilbert矩阵2—条件数 for i=3:10

h=hilb(i);

condA(i)=cond(h,2); end

disp('n cond'); for i=3:10

s=sprintf('%d %f',i,condA(i)); disp(s); end

运行后得到如下结果: n cond

3 524.056778 4 15513.738739 5 476607.250242 6 14951058.641005 7 475367356.370393 8 15257575270.772364 9 493153986466.270940 10 16025391750078.617000

结果中左边为Hilbert的阶数,右边为对应的条件数Cond2(Hn),从以上结果也可分析可知:当n的阶数稍高时, Hilbert矩阵即出现严重的病态 .

例 2.9 求

??1H2???1??21?2?? 1??3?的条件数cond2(H2),cond1(H2),cond?(H2). 解:

建立2阶Hilbert矩阵: >>H=hilb(2) H =

1.0000 0.5000 0.5000 0.3333 求其2-条件数

>>cond_2=cond(H,2) cond_2 =

35

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.3333

2.2 Matlab解线性方程组常用命令介绍

1. 求秩命令rank

在解线性方程组时,为了判断是否存在解,我们应先判断矩阵的秩.调用格式为:

c=rank(A)

2. 矩阵零空间向量null

当线性方程组Ax?0的秩r(A)?n时,方程组有无穷多个解,使用x=null(A)可得到满足Ax?0的一个解向量, 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. 其它调用格式请用help lu获得更多信息. 4. Cholsky分解chol

L=chol(A)

其中L为一个下三角矩阵,满足A?LLT. A必须为正定矩阵. 5. 条件数cond

c=cond(A,p)

A可为向量或矩阵,P取值为下列之一: 1——向量或矩阵的返回1—条件数. 2——返回向量或矩阵的2—条件数. inf——返回向量或矩阵的?—条件数.

36

'fro'——返回向量或矩阵的F—条件数. 6. 奇异值分解svd

[U,S,V]=svd(A)

将A分解为正交矩阵U,对角矩阵S和正交矩阵V的乘积,使得A=USVT. 7. 线性方程组的符号解linsolve

X=linsolve(A,b)

它等价于X=sym(A)\\sym(b).返回方程组的符号解.

3 线性方程组的迭代解法

3.1 例题解答

例 3.1 用Jacobi迭代法解以下方程:

??10x1?2x2?x3?3??2x1?10x2?x3?15 ???x1?2x2?5x3?10解:

编制迭代计算的M文件程序如下:

%Jacobi迭代法求解例3.1 % A为方程组的增广矩阵 clc;

A=[10 -2 -1 3;-2 10 -1 15;-1 -2 5 10] MAXTIME=50;%最多进行50次迭代 eps=1e-5;%迭代误差 [n,m]=size(A);

x=zeros(n,1);%迭代初值 y=zeros(n,1); k=0;

%进入迭代计算

disp('迭代过程X的值情况如下:') disp('X='); while 1

disp(x'); for i=1:1:n s=0.0; for j=1:1:n if j~=i

s=s+A(i,j)*x(j); end

y(i)=(A(i,n+1)-s)/A(i,i);

37

end end

for i=1:1:n

maxeps=max(0,abs(x(i)-y(i))); %检查是否满足迭代精度要求 end

if maxeps<=eps%小于迭代精度退出迭代 for i=1:1:n

x(i)=y(i);%将结果赋给x end return; end

for i=1:1:n%若不满足迭代精度要求继续进行迭代 x(i)=y(i); y(i)=0.0; end k=k+1;

if k>MAXTIME%超过最大迭代次数退出 error('超过最大迭代次数,退出'); return; end end

运行该程序结果如下: A =

10 -2 -1 3 -2 10 -1 15 -1 -2 5 10 迭代过程X的值情况如下:

X=

0 0 0

0.3000 1.5000 2.0000 0.8000 1.7600 2.6600 0.9180 1.9260 2.8640 0.9716 1.9700 2.9540 0.9894 1.9897 2.9823 0.9962 1.9961 2.9938 0.9986 1.9986 2.9977 0.9995 1.9995 2.9992 0.9998 1.9998 2.9997 0.9999 1.9999 2.9999 1.0000 2.0000 3.0000 1.0000 2.0000 3.0000

38

容易看出迭代计算最后结果为: x?(1,2,3)T

例 3.2 用Gauss-Seidel迭代法计算例3.1并作比较. 解:

编制求解程序Gauss_Seidel.m如下:

%Gauss_Seidel.m

%Gauss_Seidel迭代法求解例3.2 % A为方程组的增广矩阵 clc;

format long;

A=[10 -2 -1 3;-2 10 -1 15;-1 -2 5 10] [n,m]=size(A);

%最多进行50次迭代 Maxtime=50; %控制误差 Eps=10E-5; %初始迭代值 x=zeros(1,n); disp('x=');

%迭代次数小于最大迭代次数,进入迭代 for k=1:Maxtime disp(x); for i=1:n s=0.0; for j=1:n if i~=j

s=s+A(i,j)*x(j);%计算和 end end

x(i)=(A(i,n+1)-s)/A(i,i);%求出此时迭代的值 end

%因为方程的精确解为整数,所以这里将迭代结果向整数靠近的误差作为判断迭代是否停止的条件

if sum((x-floor(x)).^2)

disp('迭代结果:'); X

format short;

39

完成后直接在Matlab命令窗口中输入Gauss_Seidel回车后可得到如下结果: >>Gauss_Seidel A =

10 -2 -1 3 -2 10 -1 15 -1 -2 5 10 x=

0 0 0 0.300000000000000 2.684000000000000 0.880400000000000 2.953872000000000 0.984283200000000 2.993754176000000 0.997824185600000 2.999140939008000 0.999702144844800 2.999882238116864 0.999959128385638 2.999983845472654 0.999994394445028 2.999997784263514 0.999999231113606 2.999999696082350 0.999999894538049 2.999999958313948 0.999999985534564 2.999999994282236 0.999999998015885 2.999999999215737 0.999999999727854 2.999999999892429 0.999999999962672 2.999999999985246 0.999999999994880 2.999999999997976 0.999999999999298 2.999999999999722 0.999999999999904 2.999999999999962 0.999999999999987 2.999999999999995 0.999999999999998

40

1.560000000000000 1.944480000000000 1.992243840000000 1.998940254720000 1.999854522869760 1.999980049488814 1.999997263436271 1.999999624649072 1.999999948515845 1.999999992938308 1.999999999031401 1.999999999867145 1.999999999981778 1.999999999997501 1.999999999999657 1.999999999999953 1.999999999999994 1.999999999999999

2.999999999999999 1.000000000000000 3.000000000000000 迭代结果:

2.000000000000000

X =

1 2 3

可见对此题Gauss_Seidel法的收敛速度还是很快的.

例 3.3 取??1.45,用Gauss-Seidel迭代法和SOR法计算下列方程组的解:

?4x1?2x2?x3?0???2x1?4x2?2x3??2 ??x?2x?3x?323?1解:

Gauss-Seidel迭代法可利用上题中的程序,把输入矩阵A换掉就可以了,以下

编制求解该程序的SOR.m

%SOR法求解例3.3 %w=1.45

%方程组系数矩阵 clc;

A=[4 -2 -1;-2 4 -2 ;-1 -2 3] %方程组右端系数 b=[0,-2,3]' w=1.45;

%最大迭代次数 Maxtime=100; %精度要求 Eps=1E-5;

%以15位小数显示 format long; n=length(A); k=0;

%初始迭代值 x=ones(n,1); y=x;

disp('迭代过程:'); disp('x='); while 1 y=x; disp(x'); %计算过程 for i=1:n

41

s=b(i); for j=1:n if j~=i

s=s-A(i,j)*x(j); end end

if abs(A(i,i))<1E-10 | k>=Maxtime

error('已达最大迭代次数或矩阵系数近似为0,无法进行迭代'); return; end

s=s/A(i,i);

x(i)=(1-w)*x(i)+w*s; end

if norm(y-x,inf)

disp('最后迭代结果:'); %最后的结果 X=x'

%设为默认显示格式 format short;

为了能有可比性, Gauss?Seidel程序中的误差控制改用无穷大范数来度量(即将其改为: norm(y-x,inf)

以下为两种方法运行结果:

(1) Gauss?Seidel法运行结果如下: >>Gauss_Seidel

A =

4 -2 -1 0 -2 4 -2 -2 -1 -2 3 3 x=

1 1 1 0.750000000000000 0.375000000000000 1.500000000000000 0.562500000000000 0.531250000000000 1.541666666666667 0.651041666666667 0.596354166666667

42

1.614583333333333 0.701822916666667 1.672743055555556 0.747287326388889 1.722439236111111 0.785617404513889 1.764558015046297 0.818153663917824 1.800288447627315 0.845750031647859 1.830596170307677 0.869158662395713 1.856304498366368 0.889014832767439 1.878111387967082 0.905857679775222 1.896608915839176 0.920144495895370 1.912299302543305 0.932263178569463 1.925608553227410 0.942542758585044 1.936898023465833 0.951262333819572 1.946474230368326 0.958658646913433 1.954597174731730 0.964932513003372 1.961487400246158 0.970254271995315 1.967331981412263 0.974768413413434 1.972289602746377 0.978597499393018 1.976494867177471 0.981845492329217 1.980061950611968 0.984600577529664 1.983087701890432 0.986937557508016 1.985654272302155 0.988919882925151 1.987831346050819 0.990601375319531

0.658203125000000 0.710015190972222 0.754028320312500 0.791355839482060 0.823019239637587 0.849877416351695 0.872659665566903 0.891984533871152 0.908376705867273 0.922281240556384 0.934075655906227 0.944080178642702 0.952566438640879 0.959764843867551 0.965870836120737 0.971050197412848 0.975443551069698 0.979170179753344 0.982331264070816 0.985012629699224 0.987287077613653 0.989216360685175

43

1.989678032229960 0.992027688400078 1.991244469676705 0.993237547576686 1.992573188276692 0.994263801382521 1.993700263680578 0.995134313334948 1.994656296783491 0.995872718449754 1.995467244561000 0.996499064948561 1.996155124819374 0.997030358582234 1.996738613994614 0.997481024349056 1.997233554230909 0.997863298143645 1.997653383506066 0.998187558970155 1.998009500482125 0.998462610739587 1.998311573987100 0.998695921302203 1.998567805530502 0.998893825204951 1.998785151980135 0.999061695678897 1.998969514445976 0.999204090526252 1.999125898499494 0.999324875867931 1.999258550078451 0.999427331111469

1.999371070767130 迭代结果:

X =

0.999514237989263 1.999466515581885

(2)SOR法(??1.45)运行结果如下: >>SOR

0.990852860315019 0.992241008626696 0.993418494829607 0.994417288507763 0.995264507616623 0.995983154754781 0.996592741700804 0.997109819171835 0.997548426187277 0.997920471238110 0.998236055610856 0.998503747644651 0.998730815367726 0.998923423829516 0.999086802486114 0.999225387183712 0.999342940594960

0.999442654378196

44

A =

4 -2 -1 -2 4 -2 -1 -2 3 b = 0 -2 3 迭代过程:

x=

1 1 1 0.637500000000000 1.319906250000000 0.200426953125000 1.312280505533854 0.655033522367350 1.692284842064199 0.705846820490625 1.777193200964156 0.887273028620657 1.909222168471645 0.915403031995823 1.938503269030869 0.969682405159416 1.976329684037877 0.977544672598031 1.983982498874047 0.992436751874260 1.994143211257757 0.994223554996409 1.995924740976736 0.998182969701529 1.998589611975717 0.998531880520141 1.998972439993064 0.999573680611473 1.999666110080183 0.999628130465040 1.999741500669115 0.999901875643608 1.999922019942949 0.999905729701926 1.999934903350382

45

0.012187500000000 0.371757197265625 0.534011931458842 0.773340086195768 0.858734977660893 0.936422530391512 0.962044475111776 0.983638894760734 0.990266454150192 0.995946001166684 0.997552389716733 0.999014506686967 0.999389409429154 0.999763090152167 0.999848057258249 0.999943492726323

0.999977856325231 1.999982030452898 0.999976030495353 1.999983560868664 0.999995121966767 1.999995923188666 最后迭代结果:

0.999962179037974 0.999986613620394 0.999990568926511

X =

0.999993879742566 0.999996851108214 1.999995832511947

由以上两种方法的迭代结果可知,SOR法有明显加快收敛的作用.SOR法是广泛使用的求解大型线性方程的迭代法,只要控制好松弛因子,常可获得较好的加速速度,而目前松弛因子的选取尚无有效的方法,只能通过经验或反复的试验才能获得,这也是SOR的一个不足之处.

例3.4 就例3.1中的系数阵A1和例3.3中的系数阵A2:

?10?2?1??4?2?1?????A1???210?1?,A2???24?2?

??1?15???1?23?????讨论Jacobi迭代法和Gauss-Seidel迭代法的收敛性.

解:

先讨论它对Jacobi法的收敛性,再求解它们对的收敛性,求解程序如下:

%计算例3.4

%example3_4.m clear all; clc;

A1=[10 -2 -1;-2 10 -1;-1 -2 5] A2=[4 -2 -1;-2 4 -2;-1 -2 3] %对A1求解

%构造A1的Jocobi迭代矩阵BJ n=max(size(A1)); D=zeros(n); for i=1:n

D(i,i)=A1(i,i); end

BJ_1=-1.0*inv(D)*(A1-D)

46

%利用定理3.2判断A1对Jacobi迭代法的收敛性, %为此求其无穷大范数 BJ1_inf=norm(BJ_1,inf) %对A2求解

%构造A2的Jocobi迭代矩阵BJ n=max(size(A2)); D=zeros(n); for i=1:n

D(i,i)=A2(i,i); end

BJ_2=-1.0*inv(D)*(A2-D)

%利用定理3.2判断A2对Jacobi迭代法的收敛性, %求无穷大范数

BJ2_inf=norm(BJ_2,inf)

%以下判断对Gauss_Seidel迭代法的收敛性 %对A1

D=zeros(n); L=D; U=D;

for i=1:n

for j=1:n if i==j

D(i,i)=A1(i,i); elseif i>j

L(i,j)=A1(i,j); elseif i

U(i,j)=A1(i,j); end

end end

BS1=-1.0*inv(D+L)*U BS1_inf=norm(BS1,inf) %对A2 D=zeros(n); L=D;

47

U=D;

for i=1:n

for j=1:n if i==j

D(i,i)=A2(i,i); elseif i>j

L(i,j)=A2(i,j); elseif i

U(i,j)=A2(i,j); end

end end

BS2=-1.0*inv(D+L)*U BS2_inf=norm(BS2,inf)

程序运行结果如下: A1 =

10 -2 -1 -2 10 -1 -1 -2 5 A2 =

4 -2 -1 -2 4 -2 -1 -2 3 BJ_1 =

0 0.2000 0.1000 0.2000 0 0.1000 0.2000 0.4000 0 BJ1_inf = 0.6000 BJ_2 =

0 0.5000 0.2500 0.5000 0 0.5000 0.3333 0.6667 0 BJ2_inf = 1 BS1 =

0 0.2000 0.1000 0 0.0400 0.1200 0 0.0560 0.0680 BS1_inf = 0.3000 BS2 =

0 0.5000 0.2500

48

0 0.2500 0.6250 0 0.3333 0.5000 BS2_inf = 0.8750

Jacobi从运行结果分析可知: BJ(A1)??0.6?1,Bs(A1对1)??0.3?1,因此A迭代法和Gauss?Seidel迭代法都收剑. Bs(A2)??0.875?1,A2对Gauss?Seidel迭代法收敛,而BJ(A2)??1,故A2无法用定理3.2判断其收敛性. 为判断A2对Jacobi的收敛性,改为求其谱半径的方法. >>R_A2=max(abs(eig(BJ_2)))

R_A2 = 0.9207

其谱半径?(BJ(A2))?0.9207?1,可见A2对Jacobi迭代法还是收敛的.

例 3.5 就例3.1中系数矩阵A用从矩阵中判断收敛性的方法判断它对迭代法的收敛性,A矩阵为:

?10?2?1???A???210?1?

??1?25???解:

对该题我们可编制一个判断矩阵是否对角占优的函数IsDiagominant(A),程序编制如下：

function x=IsDiagominant(A)

%判断输入的矩阵A是否具有严格对角优势或A不可约且具有对角优势 %返回值为1,表示A是具有严格对角优势或A不可约且具有对角优势 %若返回值为0则没有上述性质 [m,n]=size(A); if m~=n

error('A不是方阵'); return; end

%先判断是否为严格对角占优

49

B=abs(A); for i=1:n b=B(i,i); temp=B(i,:); temp(i)=[]; if b>sum(temp) T(i)=1; else

T(i)=0; end end

if sum(T)==sum(ones(n,1))

Fg1=1;%输入方阵为严格对角占优矩阵 else

Fg1=0;%输入方阵产不是严格对角占优矩阵 end

%以下判断它是否为不可约且具有对角优势 %先判断是否为不可约矩阵 D=A;

D=rref(D); d=D(n-1,:); if sum(d)~=0

flag=1;%A为不可约矩阵 else

flag=0;%A为可约矩阵 end

for i=1:n b=B(i,i); temp=B(i,:); temp(i)=[];

if b>=sum(temp) T1(i)=1; else

T1(i)=0; end end

Fg2=0;%先假定A为可约矩阵且不具有对角优势

50

本文来源：https://www.bwwdw.com/article/wyw5.html