《网络数值线性代数.pdf》由会员分享，可在线阅读，更多相关《网络数值线性代数.pdf（4页珍藏版）》请在三一文库上搜索。
1、实验四实验四高斯消去法高斯消去法 4.1实验目的实验目的 1掌握高斯消去法的基本思路和迭代步骤； 2培养编程与上机调试能力。 4.2算法描述算法描述 4.2.1高斯消去高斯消去法基本思路。法基本思路。 设有方程组Axb，设A是可逆矩阵。高斯消去法的基本思想就是僵局真的初等行变 换作用于方程组的增广矩阵BA b，将其中的A变换成一个上三角矩阵，然后求解这个 三角形方程组。 4.2.2列主元高斯消去列主元高斯消去法法计算步骤计算步骤 将方程组用增广矩阵 (1) ij nn BA ba 表示。 步骤 1：消元过程，对1,2,1kn （1）选主元，找,1, k ik kn 使得 , max k iki
2、k k i n aa （2）如果 , 0 k ik a，则矩阵A奇异，程序结束；否则执行（3） 。 （3）如果 k ik， 则交换第k行与第 k i行对应元素位置， k kji j aa，,1jkn。 （4）消元，对,ikn，计算/, ikikkk laa对1,1jkn，计算 . ijijikkj aal a 步骤 2：回代过程： （1）若0, nn a则矩阵奇异，程序结束；否则执行（2） 。 （2） ,1/ ; nn nnn xaa 对1,2,1in ，计算 ,1 1 / n ii nijjii j i xaa xa 4.3实验内容实验内容 解方程组 123 123 123 0.1012.3
3、043.5551.183 1.3473.7124.6232.137 2.8351.0725.6433.035 xxx xxx xxx function X=gaussian(A,B,n) C=A B; for i=2:n r,c=max(abs(C(i-1:n,i-1); VEC=C(i-1,:);C(i-1,:)=C(c+i-2,:);C(c+i-2,:)=VEC; if C(i-1,i-1)=0 return end for k=(i-1):(n-1) C(k+1,:)=C(i-1,:)*-C(k+1,i-1)/C(i-1,i-1)+C(k+1,:); end end D=C; B2=C(
4、:,n+1); D(:,n+1)=; X=zeros(1,n); for i=1:n X(n+1-i)=(B2(n+1-i)-D(n+1-i,:)*X)/C(n+1-i,n+1-i); end = 运行: A=10,-1,2,0;-1,11,-1,3;2,-1,10,-1;0,3,-1,8; B=6,25,-11,15; gaussian(A,B,4) ans = 1.00002.0000-1.00001.0000 14 回答者：雪 function RA,RB,n,X=gaus(A,b) B=A b; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-
5、RA; if zhica0, disp(请注意：因为 RA=RB，所以此方程组无解.) return end if RA=RB if RA=n disp(请注意：因为 RA=RB=n，所以此方程组有唯一解.) X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1 for k=p+1:n m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1); end end b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(
6、A(q,q+1:n)*X(q+1:n)/A(q,q); end else disp(请注意：因为 RA=RB A=1 -1 1 -3; 0 -1 -1 1;2 -2 -4 6;1 -2 -4 1; b=1;0; -1;-1; RA,RB,n,X =gaus (A,b) 运行后输出结果 请注意：因为 RA=RB=n，所以此方程组有唯一解. RA = 4 RB = 4 n = 4 LULU 分解法分解法 .4.2.4.2直接直接 LULU 分解法及其分解法及其 MATLABMATLAB 程序程序 将矩阵将矩阵A进行直接进行直接 LU 分解分解的的 MATLABMATLAB 程序程序 functio
7、n hl=zhjLU(A) n n =size(A); RA=rank(A); if RA=n disp(请注意：因为 A 的 n 阶行列式 hl 等于零，所以 A 不能进行 LU 分 解.A 的秩 RA 如下:), RA,hl=det(A); return end if RA=n for p=1:n h(p)=det(A(1:p, 1:p); end X = 0 -0.5000 0.5000 0 hl=h(1:n); for i=1:n if h(1,i)=0 disp(请注意： 因为A 的 r阶主子式等于零， 所以A 不能进行 LU分解.A 的秩 RA 和各阶顺序主子式值 hl 依次如下：
8、), hl;RA return end end if h(1,i)=0 disp(请注意：因为 A 的各阶主子式都不等于零，所以 A 能进行 LU 分解.A 的秩 RA 和各阶顺序主子式值 hl 依次如下：) for j=1:n U(1,j)=A(1,j); end for k=2:n for i=2:n for j=2:n L(1,1)=1;L(i,i)=1; if ij L(1,1)=1;L(2,1)=A(2,1)/U(1,1); L(i,1)=A(i,1)/U(1,1); L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k)/U(k,k); else U(k,j)=
9、A(k,j)-L(k,1:k-1)*U(1:k-1,j); end end end end hl;RA,U,L end end 例例 3.4.33.4.3用矩阵进行直接 LU 分解的 MATLAB 程序分解矩阵 3010 3421 1010 0201 A . 解解 在 MATLAB 工作窗口输入程序 A=1 0 2 0;0 1 0 1;1 2 4 3;0 1 0 3; hl=zhjLU(A) 运行后输出结果 请注意：因为 A 的各阶主子式都不等于零，所以 A 能进行 LU 分解.A 的秩 RA 和各阶顺序主子式值 hl 依次如下： RA = 4 U = 1020 0101 0021 0002 L = 1000 0100 1210 0101 hl = 1124