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