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

130213037+左芝+矩阵计算上机作业七

作者:130213037左芝   发布时间:2015-06-23 17:01:01   浏览次数:87

 %列选主元

function y=gauss(A,b)
[m,n]=size(A);
p=[];
P=eye(n);
tol=1e-16;
if m==n
   for k=1:n-1
       [r,c]=find(A==max(max(abs(A(k:n,k:n)))));
       p(k)=r(1);
       A([k,p(k)],:)=A([p(k),k],:);
       P([k,p(k)],:)=P([p(k),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
           L='A(k,k)=0!';
           break;
       end
   end
else
   L='m~=n';
end
L=eye(n)+tril(A,-1);
U=triu(A,0);
b=P*b;
for j=1:(n-1)
    b(j)=b(j)/L(j,j);
    b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);
end
b(n)=b(n)/L(n,n);
y=b;
for j=n:-1:2
    y(j)=y(j)/U(j,j);
    y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);
end
y(1)=y(1)/U(1,1);

%矩阵

function [A,b,x]=matrix(n)

A=ones(n,n);

A1=(-2)*tril(A,-1);

A2=(-1)*triu(A,1);

A=A+A1+A2;

A(:,n)=1;

x=rand(n,1);

b=A*x;

 

%验证

for n=5:30

    [A,b,x]=matrix(n);

    y=gauss(A,b);

    error=norm(x-y,2)

end

 

%运行结果

error =

 

  2.2204e-016

 

 

error =

 

  1.5701e-016

 

 

error =

 

  3.2311e-015

 

 

error =

 

  5.3325e-015

 

 

error =

 

  1.3521e-015

 

 

error =

 

  4.4717e-015

 

 

error =

 

  7.7664e-015

 

 

error =

 

  2.6142e-014

 

 

error =

 

  4.3588e-014

 

 

error =

 

  2.4976e-013

 

 

error =

 

  2.8337e-013

 

 

error =

 

  8.0581e-013

 

 

error =

 

  4.1528e-013

 

 

error =

 

  5.7491e-012

 

 

error =

 

  1.3495e-012

 

 

error =

 

  2.0822e-011

 

 

error =

 

  7.5540e-012

 

 

error =

 

  1.0360e-010

 

 

error =

 

  4.0775e-011

 

 

error =

 

  1.4145e-010

 

 

error =

 

  3.6417e-010

 

 

error =

 

  1.0592e-009

 

 

error =

 

  1.0323e-009







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

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

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