矩阵分解与线性方程组求解

一、 用列主元素高斯削去法求解下述线性方程组：

?x1?13x2?2x3?34x4?13?2x?6x?7x?10x??22?1234 ??10x?x?5x?9x?141234????3x1?5x2?15x4??36程序：

function x=gaussa(a)

m=size(a); n=m(1); x=zeros(n,1); for k=1:n-1

[c,i]=max(abs(a(k:n,k))); q=i+k-1; if q~=k

d=a(q,:);a(q,:)=a(k,:);a(k,:)=d end

for i=k+1:n

a(i,:)=a(i,:)-a(k,:)*a(i,k)/a(k,k) end end

for j=n:-1:1

x(j)=(a(j,n+1)-a(j,j+1:n)*x(j+1:n))/a(j,j) end

执行过程：

>> a=[1 13 -2 -34 13;2 6 -7 -10 -22;-10 -1 5 9 14; -3 -5 0 15 -36] a =

-10 -1 5 9 14 2 6 -7 -10 -22 1 13 -2 -34 13 -3 -5 0 15 -36 >> gaussa(a) a =

-10.0000 -1.0000 5.0000 9.0000 14.0000 0 5.8000 -6.0000 -8.2000 -19.2000 1.0000 13.0000 -2.0000 -34.0000 13.0000 -3.0000 -5.0000 0 15.0000 -36.0000 a =

-10.0000 -1.0000 5.0000 9.0000 14.0000 0 5.8000 -6.0000 -8.2000 -19.2000 0 12.9000 -1.5000 -33.1000 14.4000 -3.0000 -5.0000 0 15.0000 -36.0000 a =

-10.0000 -1.0000 5.0000 9.0000 14.0000 0 5.8000 -6.0000 -8.2000 -19.2000 0 12.9000 -1.5000 -33.1000 14.4000 0 -4.7000 -1.5000 12.3000 -40.2000

a =

-10.0000 -1.0000 5.0000 9.0000 14.0000 0 12.9000 -1.5000 -33.1000 14.4000 0 5.8000 -6.0000 -8.2000 -19.2000 0 -4.7000 -1.5000 12.3000 -40.2000 a =

-10.0000 -1.0000 5.0000 9.0000 14.0000 0 12.9000 -1.5000 -33.1000 14.4000 0 0.0000 -5.3256 6.6822 -25.6744 0 a =

-10.0000 0 0 0 a =

-10.0000 0 0 0 x =

0 0 0 10.7786 x =

0 0 18.3452 10.7786 x =

0 30.9062 18.3452 10.7786 x =

14.3827 30.9062 18.3452 10.7786 ans =

14.3827 30.9062 18.3452

-4.7000 -1.0000 12.9000 0.0000 0 -1.0000 12.9000 0.0000 -0.0000 -1.5000 12.3000 -40.2000 5.0000 9.0000 14.0000 -1.5000 -33.1000 14.4000 -5.3256 6.6822 -25.6744 -2.0465 0.2403 -34.9535 5.0000 9.0000 14.0000 -1.5000 -33.1000 14.4000 -5.3256 6.6822 -25.6744 0 -2.3275 -25.0873 10.7786

二、用矩阵A的杜丽特尔三角分解A=LU，求解方程组：

?157010??x1??????618159??x2??010287??x?=

3?????50635??x????4??8????6??4? ???2???程序：

function x=lua(a,b)

m=size(a);n=m(1);x=zeros(n,1);y=zeros(n,1); for k=1:n

a(k,k:n) = a(k,k:n)-a(k,1:k-1)*a(1:k-1,k:n);

a(k+1:n,k)=(a(k+1:n,k)-a(k+1:n,1:k-1)*a(1:k-1,k))/a(k,k); end

U = triu(a,0)

L =eye(n)+ tril(a,-1) for i=1:n

y(i) = b(i)-L(i,1:i-1)*y(1:i-1); end

for i=n:-1:1

x(i)=(y(i)-U(i,i+1:n)*x(i+1:n))/U(i,i); end 执行过程：

>> a=[15 7 0 10;6 18 15 9;0 10 28 7;5 0 6 35] a =

15 7 0 10 6 18 15 9 0 10 28 7 5 0 6 35 >> b=[8;6;4;2] b = 8 6 4 2 >> lua(a,b) U =

15.0000 7.0000 0 10.0000 0 15.2000 15.0000 5.0000 0 0 18.1316 3.7105 0 0 0 30.7351 L =

1.0000 0 0 0 0.4000 1.0000 0 0

0 0.6579 1.0000 0 0.3333 -0.1535 0.4579 1.0000 ans =

0.5264 0.0718 0.1272 -0.0399

三、用乔列斯基分解计算下述线性方程组：

00??x1??5??4?10??????0??x2??8???14?10?0?14?10? ?x?=?16? ???3???0?14?1??x4??24??0?????000?14????x5??36?程序：

function x=choleskey(a,b)

m=size(a);n=m(1);x=zeros(n,1);y=zeros(n,1); G=zeros(m); for j=1:n

G(j,j)=sqrt(a(j,j)-G(j,1:j-1)*G(j,1:j-1)')

G(j+1:n,j)=(a(j+1:n,j)-G(j+1:n,1:j-1)*G(j,1:j-1)')/G(j,j) end

for i=1:n

y(i)=(b(i)-G(i,1:i-1)*y(1:i-1))/G(i,i); end

for i=n:-1:1

x(i)=(y(i)-G(i+1:n,i)'*x(i+1:n))/G(i,i); end 运行过程：

>> a=[4 -1 0 0 0;-1 4 -1 0 0;0 -1 4 -1 0;0 0 -1 4 -1;0 0 0 -1 4] a =

4 -1 0 0 0 -1 4 -1 0 0 0 -1 4 -1 0 0 0 -1 4 -1 0 0 0 -1 4 >> b=[5 ;8;16;24;36] b = 5 8 16 24 36

>> choleskey(a,b) G =

2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 G =

2.0000 0 0 0 0 -0.5000 0 0 0 G =

2.0000 -0.5000 0 0 0 G =

2.0000 -0.5000 0 0 0 G =

2.0000 -0.5000 0 0 0 G =

2.0000 -0.5000 0 0 0 G =

2.0000 -0.5000 0 0 0 G =

0 0 0 0 0 1.9365 0 0 0 0 1.9365 -0.5164 0 0 0 1.9365 -0.5164 0 0 0 1.9365 -0.5164 0 0 0 1.9365 -0.5164 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.9322 0 0 0 0 1.9322 -0.5175 0 0 0 1.9322 -0.5175 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 …… 此处隐藏：1975字，全部文档内容请下载后查看。喜欢就下载吧 ……