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

130213022_苏楠_作业7

作者:苏楠   发布时间:2015-06-23 16:31:39   浏览次数:100

1.估计计算解x的精度

程序代码:

列选主元Gauss消去法:

function [L,U,P]=gauss4(A)

n=size(A,1);

P=eye(n);

for k=1:n-1

 [r,m]=max(abs(A(k:n,k)));

 m=m+k-1;

 A([m,k],:)=A([k,m],:);

 P([m,k],:)=P([k,m],:);

 if A(k,k)~=0

    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

      stop;

    end

end

L=eye(n,n)+tril(A,-1);

U=triu(A);

运行结果:

P71

clc

format long

for n=5:30

%矩阵A的表示

A=tril(-1*ones(n),0)+2*eye(n);

A(1:n-1,n)=1;

[L,U,P]=gauss4(A);

n=size(A,1);

%任取n1列的x并求出b

x=rand(n,1);

b=A*x;

%b的无穷范数

normb=norm(b,inf);

x1=Ab;

%前代法Ly=b     

b=P*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 j=n-1:-1:1

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

    end

   x=b;

%A的转置无穷范数

norminvA=norm(inv(A),inf);

%Ax=b得到的计算解为x1

%r表示残差

r=A*x1-A*x;

normr=norm(r,inf);

normA=norm(A,inf);

%估计计算解精度

p1=norminvA*normA*normr/normb

%相对误差

error=norm((x1-x),inf)/norm(x1,inf)

%精度与相对误差的比较

err=abs(p1-error)

end

上机运行结果:

p1 =

    2.069150587705709e-015

error =

    7.074158704009518e-016

err =

    1.361734717304757e-015

p1 =

    4.283155665790917e-015

error =

    9.111116697715741e-016

err =

    3.372043996019344e-015

p1 =

    9.970535538145916e-016

error =

    1.171567088203780e-016

err =

    8.798968449942136e-016

p1 =

    1.530691775978788e-014

error =

    6.416560461151808e-015

err =

    8.890357298636070e-015

p1 =

    7.483560838186931e-015

error =

    2.054474188443684e-015

err =

    5.429086649743247e-015

p1 =

    5.926190376815989e-015

error =

    3.905365362618875e-015

err =

    2.020825014197113e-015

p1 =

    1.675342256297959e-014

error =

    4.121515135948625e-015

err =

    1.263190742703097e-014

p1 =

    1.194504638925583e-013

error =

    2.848156409099929e-014

err =

    9.096889980155898e-014

p1 =

    2.724262973554541e-013

error =

    1.151055648876753e-013

err =

    1.573207324677789e-013

p1 =

    1.854757620429866e-013

error =

    6.070295750209746e-014

err =

    1.247728045408891e-013

p1 =

    3.346845507069762e-012

error =

    4.853825906028652e-013

err =

    2.861462916466897e-012

p1 =

    7.008333894664603e-013

error =

    2.487727552827010e-013

err =

    4.520606341837594e-013

p1 =

    1.410033225864273e-011

error =

    3.737343233559194e-012

err =

    1.036298902508353e-011

p1 =

    5.818834974678200e-012

error =

    1.990849176060063e-012

err =

    3.827985798618137e-012

p1 =

    8.884517216339096e-011

error =

    2.963082749843586e-011

err =

    5.921434466495511e-011

p1 =

    3.422563790102820e-011

error =

    7.627922788400151e-012

err =

    2.659771511262805e-011

p1 =

    3.438940678012612e-011

error =

    1.516613570686042e-011

err =

    1.922327107326569e-011

p1 =

    5.927115030422014e-014

error =

    1.695258435907497e-014

err =

    4.231856594514516e-014

p1 =

    2.312049990594028e-010

error =

    6.282724539552893e-011

err =

    1.683777536638739e-010

p1 =

    1.136364437224556e-009

error =

    4.721722120234855e-010

err =

    6.641922252010702e-010

p1 =

    1.189771582851931e-009

error =

    4.951177798942828e-010

err =

    6.946538029576485e-010

p1 =

    1.350757883321099e-008

error =

    3.944095021094752e-009

err =

    9.563483812116237e-009

p1 =

    1.975082728132751e-009

error =

    4.886634507070422e-010

err =

    1.486419277425709e-009

p1 =

    3.854294710725998e-008

error =

    7.666919632831517e-009

err =

    3.087602747442846e-008

p1 =

    4.198851487065598e-008

error =

    1.529845497983454e-008

err =

    2.669005989082144e-008

p1 =

    1.742092256562085e-008

error =

    3.768168294749838e-009

err =

    1.365275427087101e-008

2.利用2.5.2迭代改进算法计算x的精度

程序代码:

Gauss消去法:

function [L,U]=gauss(A)

n=size(A,1);

for k=1:n-1

    if A(k,k)~=0

       tol=10e-16;

       for k=1:n-1

       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

          stop;

       end

       end

    end

end

L=eye(n)+tril(A,-1);

U=triu(A,0);

//

运行结果:

P72  2.5.2迭代改进

clc

format long

for n=5:30

 %矩阵A的表示

A=tril(-1*ones(n),0)+2*eye(n);

A(1:n-1,n)=1;

[L,U,P]=gauss4(A);

n=size(A,1);

%任取n1列的x并求出b

x1=rand(n,1);

d=A*x1;

b=d;

%b的无穷范数

normb=norm(d,inf);

x1=Ab;

%前代法Ly=b

b=P*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 j=n-1:-1:1

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

    end

   x=b;

%A的转置无穷范数

norminvA=norm(inv(A),inf);

%Ax=b得到的计算解为x1

%r表示残差

r=A*x1-A*x;

normr=norm(r,inf);

normA=norm(A,inf);

%估计计算解精度

p1=norminvA*normA*normr/normb;

%相对误差

error=norm((x1-x),inf)/norm(x1,inf);

%精度与相对误差的比较

err=abs(p1-error);

while error>1e-16

r1=d-A*x;

[L,U]=gauss(A);

n=size(A,1);

%前代法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 j=n-1:-1:1

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

    end

   z=b;

  x1=x+z;

   x=x1;

%估计计算解精度

normr1=norm(r1,inf);

p2=norminvA*normA*normr1/normb;

%相对误差

error=norm((x1-x),inf)/norm(x1,inf);

%精度与相对误差的比较

er=abs(p2-error)

end

end

上机运行结果:

er =

    5.157203468659786e-017

er =

    1.272496802478046e-014

er =

    2.155269647583724e-014

er =

    1.871617510843322e-015

er =

    5.951148823082541e-014

er =

    1.902524521511317e-013

er =

    3.163071505155942e-013

er =

    1.710248286968516e-013

er =

    4.973030961919993e-013

er =

    9.533093019739630e-013

er =

    3.867972841067116e-012

er =

    1.740282290643754e-012

er =

    9.587851737573304e-012

er =

    3.611249303813100e-011

er =

    1.384995088127064e-011

er =

    1.527385147351999e-010

er =

    1.348099132693654e-010

er =

    8.610359087827285e-011

er =

    4.103438333102761e-010

er =

    8.437357824666399e-010

er =

    6.300959147279881e-010

er =

    1.185281282445974e-009

er =

    1.194010099955666e-008

er =

    4.742954132219643e-008

er =

    7.724249723540388e-010







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

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

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