现在时间是:
当前位置:首 页 >> 矩阵计算>> 教学区>> 文章列表

130213004 葛嫚嫚 全选主元与列选主元消去法

作者:葛嫚嫚   发布时间:2015-06-23 16:25:34   浏览次数:117

function x=jiex(A,b)

format long

A=[0.001 1.00;1.00 2.00];

b=[1.00;3.00];

[L,U]=gass(A);

y=qiandaiy(L,b);

x=huidaix(U,y);

 

function [L,U]=gass(A)

n=size(A,1)

tol=10^(-16);

for k=1:n-1

  if abs(A(k,k))>tol

      A(k+1:n,k)=A(k+1:n,k)/A(k,k);

      A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);

  end

end

  L=eye(n)+tril(A,-1);

  U=triu(A,0);

 

function y=qiandaiy(L,b)

n=size(L,1);

y=zeros(n,1);

y(1)=b(1)/L(1,1);

for j=2:n

  y(j)=(b(j)-L(j,1:j-1)*y(1:j-1))/L(j,j);

end

 

function x=huidaix(U,y)

n=size(U,1);

x=zeros(n,1);

x(n)=y(n)/U(n,n);

for j=n-1:-1:1

x(j)=(y(j)-U(j,j+1:n)*x(j+1:n))/U(j,j);

end


结果:ans =

  1.00200400801598

  0.99899799599198

精确值:x=inv(A)*b

inv(A)*b

ans =

  1.00200400801603

  0.99899799599198


全选主元

(2.1)

function x=jiequanx(A,b)

A=[0.001 1.00;1.00 2.00];

b=[1.00;3.00];

[L,U,p,q]=quan(A);

for k=1:length(p)

   b([k,p(k)])=b([p(k),k]);  %交换k,p(k)列

end             

y=qiandaiy(L,b); 

x=huidaix(U,y); 

for k=length(q):-1:1  

   x([k,q(k)])=x([q(k),k]);  %交换k,q(k)列

end


function [L,U,p,q]=quan(A,b)

[m,n]=size(A);

tol=10e-16;

p=[];q=[];

if m==n

for k=1:n-1

[r,c]=find(A==max(max(A(k:n,k:n))));  %找到A中最大数的位置

p(k)=r

q(k)=c

A([k,p(k)],:)=A([p(k),k],:);   %交换k行和p(k)行

A(:,[k,q(k)])= A(:,[q(k),k]);  %交换k列和q(k)列

 if abs(A(k,k))>tol

  A(k+1:n,k)=A(k+1:n,k)/A(k,k);

  A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);

 else

     'A(k,k)=0!';

     break;

 end

end

else 'm~=n';

end

L=eye(n)+tril(A,-1);

U=triu(A,0);


function y=qiandaiy(L,b)

n=size(L,1);

y=zeros(n,1);

y(1)=b(1)/L(1,1);

for j=2:n

  y(j)=(b(j)-L(j,1:j-1)*y(1:j-1))/L(j,j);

end

 

function x=huidaix(U,y)

n=size(U,1);

x=zeros(n,1);

x(n)=y(n)/U(n,n);

for j=n-1:-1:1

x(j)=(y(j)-U(j,j+1:n)*x(j+1:n))/U(j,j);

end

 

结果:ans =

 1.002004008015978

 0.998997995991984

精确值:x=inv(A)*b

ans =

  1.00200400801603

  0.99899799599198


列选主元

(2.2) function x=jieliex(A,b)

A=[0.001 1.00;1.00 2.00];

b=[1.00;3.00];

[L,U,p]=lie(A);

for k=1:length(p)

b([k,p(k)])=b([p(k),k]);   %交换k,p(k)列

end

y=qiandaiy(L,b);

x=huidaix(U,y);


function [L,U,p]=lie(A)

n=size(A,1);

tol=1e-16;

P=[];

for k=1:n-1

[r,c]=find(A==max(max(abs(A(k:n,k:n)))));  %找到A中最大数的位置

p(k)=r;    

A([p(k),k],:)=A([k,p(k)],:);    %交换k行和p(k)行

if abs(A(k,k))>tol

A(k+1:n,k)=A(k+1:n,k)/A(k,k);

A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);

end

end

L=eye(n)+tril(A,-1);

U=triu(A,0);


function y=qiandaiy(L,b)

n=size(L,1);

y=zeros(n,1);

y(1)=b(1)/L(1,1);

for j=2:n

y(j)=(b(j)-L(j,1:j-1)*y(1:j-1))/L(j,j);

end

function x=huidaix(U,y)

n=size(U,1);

x=zeros(n,1);

x(n)=y(n)/U(n,n);

for j=n-1:-1:1

x(j)=(y(j)-U(j,j+1:n)*x(j+1:n))/U(j,j);

end

结果:ans =


1.00200400801603

0.99899799599198

精确值:x=inv(A)*b

ans =

1.00200400801603

0.99899799599198

 







上一篇:没有了    下一篇:没有了

Copyright ©2019    计算数学达人 All Right Reserved.

技术支持:自助建站 | 领地网站建设 |短信接口 |燕窝 版权所有 © 2005-2019 lingw.net.粤ICP备16125321号 -5