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

130213030严玉婷作业7

作者:严玉婷   发布时间:2015-06-23 16:31:40   浏览次数:153

%求矩阵A,随机选取x1,求b=A*x1
function [A,b]=qiu_Ab(n)
B=-1*ones(n,n);
C=zeros(n,n);
for i=1:n-1
    C(i,n)=1;
end
A=C+tril(B,-1)+eye(n); %表示出A
x1=rand(n,1);
b=A*x1;    %计算出b

 


%用列选主元的方法求解x
function x=jieliex(A,b)
[L,U,p]=lie(A);
for k=1:length(p)
b([k,p(k)])=b([p(k),k]);
end
y=qiandaiy(L,b);
x=huidaix(U,y);

 


%列选主元 LU 分解PA=LU
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)))));
%p(k)=r;
%A([p(k),k],:)=A([k,p(k)],:);
[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);
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

 


%计算解的精度,相对误差,比较两者
for n=5:1:30
[A,b]=qiu_Ab(n);
x=jieliex(A,b);
r=b-A*x;
jingdu=cond(A,inf)*norm(r,inf)/norm(b,inf)%解的精度
xiangdui_err=norm(r,inf)/norm(b,inf)%相对误差
bijiao=norm(jingdu-xingdui_err)
end

 

结果:
jingdu =

  4.7935e-016


xiangdui_err =

  9.5870e-017


bijiao =

  2.3184e-009


jingdu =

  3.6519e-015


xiangdui_err =

  6.0864e-016


bijiao =

  2.3184e-009


jingdu =

  1.5472e-015


xiangdui_err =

  2.2103e-016


bijiao =

  2.3184e-009


jingdu =

  3.6644e-014


xiangdui_err =

  4.5806e-015


bijiao =

  2.3184e-009


jingdu =

  2.4627e-014


xiangdui_err =

  2.7364e-015


bijiao =

  2.3184e-009


jingdu =

  1.3074e-013


xiangdui_err =

  1.3074e-014


bijiao =

  2.3183e-009


jingdu =

  2.9106e-014


xiangdui_err =

  2.6460e-015


bijiao =

  2.3184e-009


jingdu =

  8.5027e-014


xiangdui_err =

  7.0856e-015


bijiao =

  2.3183e-009


jingdu =

  2.2656e-013


xiangdui_err =

  1.7428e-014


bijiao =

  2.3182e-009


jingdu =

  3.0023e-013


xiangdui_err =

  2.1445e-014


bijiao =

  2.3181e-009


jingdu =

  1.0058e-012


xiangdui_err =

  6.7053e-014


bijiao =

  2.3174e-009


jingdu =

  3.2843e-012


xiangdui_err =

  2.0527e-013


bijiao =

  2.3151e-009


jingdu =

  1.7903e-013


xiangdui_err =

  1.0531e-014


bijiao =

  2.3183e-009


jingdu =

  7.4203e-012


xiangdui_err =

  4.1224e-013


bijiao =

  2.3110e-009


jingdu =

  4.3398e-012


xiangdui_err =

  2.2841e-013


bijiao =

  2.3141e-009


jingdu =

  2.8784e-011


xiangdui_err =

  1.4392e-012


bijiao =

  2.2896e-009


jingdu =

  1.7365e-010


xiangdui_err =

  8.2689e-012


bijiao =

  2.1448e-009


jingdu =

  1.7354e-011


xiangdui_err =

  7.8882e-013


bijiao =

  2.3011e-009


jingdu =

  2.3950e-010


xiangdui_err =

  1.0413e-011


bijiao =

  2.0789e-009


jingdu =

  7.4585e-010


xiangdui_err =

  3.1077e-011


bijiao =

  1.5726e-009


jingdu =

  4.2521e-010


xiangdui_err =

  1.7008e-011


bijiao =

  1.8932e-009


jingdu =

  2.2456e-009


xiangdui_err =

  8.6369e-011


bijiao =

  7.2843e-011


jingdu =

  1.7878e-008


xiangdui_err =

  6.6214e-010


bijiao =

  1.5559e-008


jingdu =

  8.1957e-009


xiangdui_err =

  2.9270e-010


bijiao =

  5.8773e-009


jingdu =

  9.3810e-009


xiangdui_err =

  3.2348e-010


bijiao =

  7.0626e-009


jingdu =

  5.3336e-008


xiangdui_err =

  1.7779e-009


bijiao =

  5.1017e-008

迭代改进:[A,b]=qiu_Ab(n);
x=jieliex(A,b);
x1=ones(n,1);
r=b-A*x;
[L,U]=gass(A);
y=qiandaiy(L,b);
z=huidaiz(U,y);
x1=x+z;
if(norm((x-x1),inf)/norm(x1,inf)>eps)
x=x1;
r=b-A*x;
[L,U]=gass(A);
y=qiandaiy(L,b);
z=huidaiz(U,y);
x1=x+z;
end







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

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

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