﻿ 130213022_苏楠_作业7 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 130213022_苏楠_作业7

1.估计计算解x的精度

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