matlab解线性方程组线性方程组及MATLAB应用
1matlab 解线性方程组 线性方程组及 MATLAB 应用数值实验 线性方程组与 MATLAB 应用王1．实验目的：理解矩阵的范数与条件数。 实验内容：已知矩阵111111A11111111求 A1，A2，A 和 cond2(A)。 11解：编写了一个 M 文件来求矩阵 A 的范数与条件数：test3_1.m 如下：A=[1 1 1 1;-1 1 -1 1;-1 -1 1 1;1 -1 -1 1]; norm(A,1) norm(A,2) norm(A,inf) cond(A,2)计算结果依次是： 4 2 4 1.00002．实验目的：研究高斯消去法的数值稳定性(出现小主2元) 。实验内容：设方程组 Axb，其中两个矩阵如下，分别对以上两个方程组0.310155.291(1)A111.21159.1746.786.13012 ，b11952221159.14371032.09999999999999(2)A251100185.9000000000000162 ，b25511023(1)计算矩阵的条件数，判断系数矩阵是良态的还是病态的？解： 本题编写了一个 test3_21 的 M 文件如下：A1=[0.3*1e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; A2=[10 -7 0 1;-3 2.099999999999999 6 2;5 -1 5 -1;0 1 0 2]; cond(A1) cond(A2)求得两个矩阵的条件数分别为 68.4296 和 8.9939，易知这矩阵 A1 的条件数远远大于 1，而矩阵 A2 的条件数刚大于1，故这,矩阵 A1 为病态矩阵，矩阵 A2 为良态矩阵。 (2)用列主元消去法求得 L 和 U 及解向量 x1,x2R；解：本题利用 Matlab 的列主元三角分解函数 lu()；具体求解如下： >> A1=[0.3*1e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; >> A2=[10 -7 0 1;-3 2.099999999999999 6 2;5 -1 5 -1;0 1 0 2];4>> b1=[59.17;46.78;1;2];>> b2=[8;5.0000000000001;5;1]; >> [L1,U1]=lu(A1)L1 = 0.0000 1.0000 0 0 0.4724 -0.1755 1.0000 0 1.0000 0 0 0 0.0893 0.0202 -0.1738 1.0000 U1 = 11.2000 9.0000 5.0000 2.0000 0 59.1400 3.0000 1.0000 0 0 4-2.8354 1.2307 0 0 0 1.0151 >> [L2,U2]=lu(A2)L2 =1.0000 0 0 0 -0.3000 -0.0000 1.0000 0 0.5000 1.0000 0 0 0 0.4000 -0.3333 1.0000 U2 =10.0000 -7.0000 0 1.0000 0 2.5000 5.0000 -1.5000 0 0 6.0000 2.3000 0 0 0 3.3667 >> y1=L1\b1; >> x1=U1\y1 x1 =3.8457 1.6095 -15.4761 10.4113 >> y2=L2\b2; >> x2=U2\y2 x2 =0.1337 -0.8218 0.8842 0.91091,x2R4x 用不选主元的高斯消去法求得 L 和 U 及解向量；解：编写一个 LU_Fact 的 M 文件储存不选主元的 LU 分解法然后调用求解。具体实现如下LU_Fact.m 如下所示：function [x,L,U,index]=LU_Fact(A,b)%普通 LU 分解。 A 为要分解的矩阵， b 为方程组右端项；%x 为方程组的解，L 为单位下三角阵，U 为上三角阵，index 为指示变量 [n,m]=size(A); nb=length(b); if n ~= m5error(‘矩阵 A 的行列必须相等！’);return; endif m ~= nberror(‘矩阵 A 的行数必须和 b 的行数相同！’); return; endL=eye(n);U=zeros(n);index=1;x=zeros(n,1);y=zeros(n,1); for k=1:n for j=k:n z=0;for q=1:k-1z=z+L(k,q)*U(q,j); endU(k,j)=A(k,j)-z; endif abs(U(k,k))for i=k+1:n z=0;for q=1:k-1z=z+L(i,q)*U(q,k); endL(i,k)=(A(i,k)-z)/U(k,k); end endy(1)=b(1); for k=2:n z=0;for j=1:k-1z=z+L(k,j)*y(j); endy(k)=b(k)-z; endx(n)=y(n)/U(n,n); for k=n-1:-1:1 z=0;for j=k+1:nz=z+U(k,j)*x(j); endx(k)=(y(k)-z)/U(k,k); end求解过程及结果如下： 对于方程组 1 的求解过程： >> 6[x,L,U,index]=LU_Fact(A1,b1)x = 0 0 0 0L = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1U = 0.0000 59.1400 3.0000 1.0000 0 0 0