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

130213013+李艳+矩阵计算+作业7

作者:130213013   发布时间:2015-06-23 17:25:47   浏览次数:130

 

1.Gauss.m
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.m
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.m
for n=5:30
    [A,b,x]=matrix(n);
    y=gauss(A,b);
    error=norm(x-y,2)
end
 
结果:
error =
 
 2.4825e-016
error =
 
 1.1116e-015
error =
 
 1.6120e-015
error =
 
 3.4837e-016
error =
 
 3.0860e-015
error =
 
 6.1991e-015
error =
 
 2.2847e-014
error =
 
 3.9232e-014
error =
 
 1.6496e-013
error =
 
 4.3955e-014
error =
 
 2.1440e-013
error =
 
 4.3213e-013
error =
 
 1.3134e-012
error =
 
 1.9521e-012
error =
 
 9.1956e-013
error =
 
 1.2180e-011
error =
 
 1.1582e-011
error =
 
 3.3860e-011
error =
 
 1.8725e-010
error =
 
 6.1072e-011
error =
 
 5.0874e-010
error =
 
 7.8975e-010
error =
 
 1.9822e-009
error =
 
 2.5387e-009
error =
 
 3.8170e-009
 
error =
 
 7.9459e-009
 
 
2.gaussliexuan.m
 function [L,U,q]=guassliexuan(A)
[m,n]=size(A);
p=[];
tol=1e-16;
if m==n
    for k=1:n-1
        [r,c]=max(abs(A(k:n,k)));
        q(k)=c+k-1;
       A([k,q(k)],:)=A([q(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
            'A(k,k)=0!'
            break;
        end
    end
else
    'm~=n';
end
L=eye(n)+tril(A,-1);
U=triu(A);
Test2.m
clc
error1=[];
err1=[];
k=[5:30];
for n=k
    [A,b,x]=matrix(n);
    y=gauss(A,b);%计算解
    norma=norm(inv(A),inf);
    rnorm=norm(b-A*y,inf);
    bnorm=norm(b,inf);
    anorm=norm(A,inf);
    error=(norma*anorm*rnorm)/bnorm;%计算解的精度
    error1=[error1;error];
    err=norm(abs(x-y),inf)/norm(x,inf);%真是相对误差
    err1=[err1;err];
end
error1
err1
error=abs(error1-err1)
 
 
结果:
error1 =
 
 1.0e-007 *
 
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0001
    0.0003
    0.0000
    0.0002
    0.0009
    0.0009
    0.0009
    0.0040
    0.0371
    0.0856
    0.0367
    0.0730
    0.2688
 
 
err1 =
 
 1.0e-008 *
 
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0001
    0.0001
   0.0002
    0.0005
    0.0001
    0.0008
    0.0043
    0.0023
    0.0020
    0.0117
    0.0511
    0.3175
    0.1402
    0.1343
    0.4631
 
 
error =
 
 1.0e-007 *
 
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
   0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0003
    0.0000
    0.0002
    0.0005
    0.0007
    0.0007
    0.0028
    0.0320
    0.0538
    0.0226
    0.0596
    0.2225






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

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

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