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

130213033张艳作业七

作者:130213033   发布时间:2015-06-23 16:22:48   浏览次数:97

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

format long

x=rand(n,1)

N=-1*ones(n,n);

A=tril(N);

for i=1:n

   A(i,i)=1;

end

A(:,n)=1;

b=A*x

列选主元Gauss消去法

function [L,U,q]=guass0(A)

% Gaussion Elimination with complete pivoting

[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],:); %交换k,p

        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);

 

%test_guass0

function x1=test_gauss0(A,b)

[L,U,q]=guass0(A);

n=size(A,1);

% the interchanges of b by p

for k=1:length(q)

   b([k,q(k)])=b([q(k),k]);

end

%Ly=b

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

for j=2:n

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

end

%Ux=y

b(n)=b(n)/U(n,n);

for i=n-1:-1:1

   b(i)=(b(i)-U(i,i+1:n)*b(i+1:n))/U(i,i);

end

x1=b;%计算解

clc

clear all

format long

nn=[5:30];

error=[];

terror=[];

error1=[];

terror1=[];

for n=nn

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

   normb=norm(b,inf);

   x1=test_gausslu_pp(A,b);

   

   r=A*x-A*x1;  

   normr=norm(r,inf);  

   condA=cond(A,inf);

   error=condA*(normr/normb);

   error1=[error1;error];

   terror=norm((abs(x1-x)),inf)/norm(x,inf);  

   terror1=[terror1;terror];

end

error1

terror1

err=abs(error1-terror1)

 







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

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

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