数值分析第五章

第一题： 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 =