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

陈思130213001作业七

作者:陈思   发布时间:2015-06-23 17:24:03   浏览次数:78

    

  
Gauss: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);
matrix: 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;
test: for n=5:30
    [A,b,x]=matrix(n);
    y=gauss(A,b);
    error=norm(x-y,2)
end

结果:error =

  5.1179e-016


error =

  3.0375e-016


error =

  2.3149e-016


error =

  2.9710e-015


error =

  5.5174e-015


error =

  1.7947e-015


error =

  1.3352e-014


error =

  1.4152e-014


error =

  5.2444e-014


error =

  1.9273e-013


error =

  4.5993e-013


error =

  7.4337e-013


error =

  1.4850e-012


error =

  2.0571e-012


error =

  6.4771e-012


error =

  1.0776e-012


error =

  7.8049e-012


error =

  4.6851e-011


error =

  2.4160e-011


error =

  2.5830e-011


error =

  1.2246e-010


error =

  7.5735e-010


error =

  3.2478e-009


error =

  1.5314e-009


error =

  1.9735e-009


error =

  6.1606e-009
 







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

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

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