数值分析第五章
第一题: LU分解法: 建立m文件
function h1=zhijieLU(A,b)%h1各阶主子式的行列式值 [n n]=size(A);RA=rank(A); if RA~=n
disp('请注意:因为A的n阶行列式h1等于零,所以A不能进行LU分解。A的秩RA如下:') RA,h1=det(A); return end if RA==n for p=1:n
h(p)=det(A(1:p,1:p)); end
h1=h(1:n); for i=1:n if h(1,i)==0
disp('请注意:因为A的r阶主子式等于零,所以A不能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下:') h1;RA return end end if h(1,i)~=0
disp('请注意:因为A的r阶主子式都不等于零,所以A能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下:') 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 i>j
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)=A(k,j)-L(k,1:k-1)*U(1:k-1,j); end end end end
h1;RA,U,L,X=inv(U)*inv(L)*b
end end 输入:
>> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2]; >> b=[8;5.900001;5;1]; >> h1=zhijieLU(A,b) 输出:
请注意:因为A的r阶主子式都不等于零,所以A能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下: RA = 4 U =
10.0000 -7.0000 0 1.0000 0 2.1000 6.0000 2.3000 0 0 -2.1429 -4.2381 0 -0.0000 0 12.7333 L =
1.0000 0 0 0 -0.3000 1.0000 0 0 0.5000 1.1905 1.0000 -0.0000 0.2000 1.1429 3.2000 1.0000 X =
-0.2749 -1.3298 1.2969 1.4398 h1 =
10.0000 -0.0000 -150.0001 -762.0001 列主元高斯消去法: 建立m文件
function [RA,RB,n,X]=liezhu(A,b)
B=[A b];n=length(b);RA=rank(A);RB=rank(B);zhicha=RB-RA; if zhicha>0
disp('请注意:因为RA~=RB,所以方程组无解') return
warning offMATLAB:return_outside_of_loop end if RA==RB if RA==n
disp('请注意:因为RA=RB,所以方程组有唯一解') X=zeros(n,1);C=zeros(1,n+1); for p=1:n-1
[Y,j]=max(abs(B(p:n,p)));C=B(p,:); B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;
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(A(q,q+1:n)*X(q+1:n)))/A(q,q); end else
disp('请注意:因为RA=RB 输入: >> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2]; >> b=[8;5.900001;5;1]; >> [RA,RB,n,X]=liezhu(A,b),H=det(A) 输出: 请注意:因为RA=RB,所以方程组有唯一解 RA = 4 RB = 4 n = 4 X = 0.0000 -1.0000 1.0000 1.0000 H = -762.0001 第二题: 建立列主元高斯消去法m文件(题一中已有) (1)输入: >> format compact >> A=[3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34]; >> b=[1;1;1]; >> [RA,RB,n,X]=liezhu(A,b),h=det(A),C=cond(A) 输出: 请注意:因为RA=RB,所以方程组有唯一解 RA = 3 RB =