实验五 线性代数方程组的数值解matlab代码

例： A=[1 2 3;4 5 6;7 8 9]; t=trace(A),d=diag(A) l=tril(A),u=triu(A) b=flipud(A),c=fliplr(A)
MATLAB处理稀疏矩阵: 进行大规模计算的优点
实际应用中， 实际应用中，许多大型矩阵都含有很 元素，称为稀疏矩阵。 多0元素，称为稀疏矩阵。 a=sparse(r,c,v,m,n) 在第 行、第c列 在第r行 列 输入非0元素v，矩阵共m行 列 输出a为稀 输入非0元素 ，矩阵共 行n列，输出 为稀 疏矩阵，只给出(r,c)及v 疏矩阵，只给出 及 a=sparse(aa) 将满元素矩阵 转化为稀 将满元素矩阵aa转化为稀 疏矩阵a 疏矩阵 aa=full(a) 输入稀疏矩阵 ，输出aa为满 输入稀疏矩阵a，输出 为满 矩阵（包含零元素） 矩阵（包含零元素）
线性方程组数值解法的MATLAB实现
3. 范数 条件数 特征值 n=norm (x, p) 输入x为向量或矩阵，输出为x的p-范 数,p可取1，2，‘inf’，分别为1-范数， 2-范数，∞-范数，p缺省值为2. c=cond (x，p) 输入x为矩阵, 输出为x的p-条件数，p 的含义同norm r=rcond (x) 输入x为方阵, 输出为x条件数倒数 e=eig(x) 输入x为矩阵，输出x的全部特征值
例：a=sparse(2,2:3,8,2,4), aa=full(a), a =(2,2) 8 输出 (2,3) 8 aa= 0 0 0 0 0 8 8 0
例：s=sparse([1 6 2],[2 4 6],[0.1 0.2 0.3],6 ,6) A=full(s)
例. 分别用稀疏矩阵和满矩阵求解Ax=b, 比较计算时间 分别用稀疏矩阵和满矩阵求解
其他相关的MATLAB函数
提取（产生） 1. 提取（产生）对角阵 v=diag(x) 输入向量 ，输出 是以 为对角元素的对角阵； 输入向量x，输出v是以 为对角元素的对角阵； 是以x为对角元素的对角阵 输入矩阵x，输出v是 的对角元素构成的向量 的对角元素构成的向量； 输入矩阵 ，输出 是x的对角元素构成的向量； trace(x) 返回 的迹，即对角元素的和 返回x的迹 的迹， 提取上（ 2. 提取上（下）三角阵 y=triu(x) v=tril(x) v=triu(x,1) v=tril(x,-1) flipud(x) fliplr(x) ) 输入矩阵x，输出 是 的上三角阵 的上三角阵； 输入矩阵 ，输出v是x的上三角阵； 输入矩阵x，输出v是 的下三角阵 的下三角阵； 输入矩阵 ，输出 是x的下三角阵； 同上，但对角元素为0， 同上，但对角元素为 ， 同上，但对角元素为0。 同上，但对角元素为 。 矩阵x上下翻转 矩阵 上下翻转 矩阵x 矩阵x左右翻转
例：
clear，clc A=[1 1/4 0;0 1/2 0;0 1/4 0]; x=[0 9 0.1 0]; [norm(x),norm(x,1),norm(x,inf)] [norm(A),norm(A,1),norm(A,inf)] [cond(A),cond(A,1),cond(A,inf)] rcond(A)
线性方程组数值解法的MATLAB实现
1. 求解 Ax = b 用左除： x = A\ b 。 用左除：
若A为可逆方阵，输出原方程的解x 若A为n×m矩阵（n>m），且ATA可逆，输出原 方程的最小二乘解x
线性方程组数值解法的MATLAB实现
2. 矩阵 分解 矩阵LU分解
[x,y]=lu(A) 若A可逆且顺序主子式不为零， 输出 为单位下 可逆且顺序主子式不为零， 可逆且顺序主子式不为零 输出x为单位下 三角阵L， 为上三角阵 为上三角阵U， 可逆， 为一交 三角阵 ，y为上三角阵 ，使A=LU； 若A可逆，x为一交 ； 可逆 换阵与单位下三角阵之积. 换阵与单位下三角阵之积 [x,y,p]=lu(A) 若A可逆，输出 为单位下三角阵 ， y为上三 可逆， 为单位下三角阵L， 为上三 可逆 输出x为单位下三角阵 角阵U，p为一交换阵 ，使PA=LU. 为一交换阵P， 角阵 ， 为一交换阵 u =chol(A) 对正定对称矩阵 的Cholesky分解，输出 为上 对正定对称矩阵A的 分解， 分解 输出u为上 三角阵U， 三角阵 ，使A=UTU [Q,R]=qr(A) 正交三角分解，其中，Q为正交阵，R为上三角阵 正交三角分解，其中， 为正交阵 为正交阵， 为上三角阵 [U,S,V]=svd(A) 奇异值分解，即A=USV ’ ，U，V为正交阵，S 奇异值分解， 为正交阵， ， 为正交阵 为对角阵。 为对角阵。
t1, t2相差巨大，说明用稀疏矩阵计算的优点 相差巨大， 相差巨大 说明用稀疏矩阵计算的优点 用于简单地验证两种方法结果的一致） （y=yy 用于简单地验证两种方法结果的一致）
例. 解
10x1 + 3x2 + x3 =14
并对系数矩阵 2x1 −10x2 + 3x3 =−5 作LU分解 分解
x1 + 3x2 +10x3 =14
A=[10 3 1;2 -10 3;1 3 10], b=[14 -5 14]', x=A\b, [L1,U1]=lu(A); L1,U1, A1=L1*U1, [L2,U2,P]=lu(A); L2,U2,P, A2=L2*U2, A3=inv(P)*A2 [Q,R]=qr(A) Q*R [U,S,V]=svd(A) U*S*V’
4 1 1 4 1 0 设 A= 1 4 O 0 O O 1 1 4n×n
b = [ ,2,L ] 1 n
T
n=500;b=[1:n]'; a1=sparse(1:n,1:n,4,n,n); a2=sparse(2:n,1:n-1,1,n,n); a=a1+a2+a2'; tic;x=a\b;t1=toc aa=full(a); tic;xx=aa\b;t2=toc y=sum(x) yy=sum(xx)
线性方程组数值解法的MATLAB实现
4． ilbert 矩 ： h=hilb(n) H ． 阵 输 出h为n阶Hilbert 矩 阵
L 1 1/ 2 L 1/ n 1/ 2 1/ 3 L 1/(n +1) L H= L L 1/ n 1/(n +1) L /(2n −1) 1
当n很大时Hilbert 矩阵呈病态
例：观察Hilbert矩阵的病态性 观察 矩阵的病态性
例. Hx=b, 其中 H=hilb(5), b=[1,…1]T
H=hilb(5), h=rats(H), b=ones(5,1); x=H\b; b(5)=1.1; x1=H\b; [x,x1], n1=cond(H), n2=rcond(H), x x1 1.0e+003 *Fra Baidu bibliotek0.0050 0.0680 -0.1200 -1.3800 0.6300 6.3000 -1.1200 -9.9400 0.6300 5.0400 cond(H)=4.7661e+005