文档库

最新最全的文档下载
当前位置:文档库 > 数值分析(第五版)计算实习题第五章作业

数值分析(第五版)计算实习题第五章作业

数值分析第五章

第一题:

LU分解法:

建立m文件

function h1=zhijieLU(A,b)%h1各阶主子式的行列式值

[n n]=size(A);RA=rank(A);

if RA~=n

disp('请注意:因为A的n阶行列式h1等于零,所以A不能进行LU分解。A的秩RA如下:')

RA,h1=det(A);

return

end

if RA==n

for p=1:n

h(p)=det(A(1:p,1:p));

end

h1=h(1:n);

for i=1:n

if h(1,i)==0

disp('请注意:因为A的r阶主子式等于零,所以A不能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下:')

h1;RA

return

end

end

if h(1,i)~=0

disp('请注意:因为A的r阶主子式都不等于零,所以A能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下:')

for j=1:n

U(1,j)=A(1,j);

end

for k=2:n

for i=2:n

for j=2:n

L(1,1)=1;L(i,i)=1;

if i>j

L(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1);

L(i,k)=(A(i,k)-L(i,1:k-1)*U(1:k-1,k))/U(k,k);

else

U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);

end

end

end

end

h1;RA,U,L,X=inv(U)*inv(L)*b

end

end

输入:

>> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];

>> b=[8;5.900001;5;1];

>> h1=zhijieLU(A,b)

输出:

请注意:因为A的r阶主子式都不等于零,所以A能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下:

RA =

4

U =

10.0000 -7.0000 0 1.0000

0 2.1000 6.0000 2.3000

0 0 -2.1429 -4.2381

0 -0.0000 0 12.7333

L =

1.0000 0 0 0

-0.3000 1.0000 0 0

0.5000 1.1905 1.0000 -0.0000

0.2000 1.1429 3.2000 1.0000

X =

-0.2749

-1.3298

1.2969

1.4398

h1 =

10.0000 -0.0000 -150.0001 -762.0001

列主元高斯消去法:

建立m文件

function [RA,RB,n,X]=liezhu(A,b)

B=[A b];n=length(b);RA=rank(A);RB=rank(B);zhicha=RB-RA;

if zhicha>0

disp('请注意:因为RA~=RB,所以方程组无解')

return

warning offMATLAB:return_outside_of_loop

end

if RA==RB

if RA==n

disp('请注意:因为RA=RB,所以方程组有唯一解')

X=zeros(n,1);C=zeros(1,n+1);

for p=1:n-1

[Y,j]=max(abs(B(p:n,p)));C=B(p,:);

B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;

for k=p+1:n

m=B(k,p)/B(p,p);

B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);

end

end

b=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);

for q=n-1:-1:1

X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end

else

disp('请注意:因为RA=RB

end

输入:

>> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];

>> b=[8;5.900001;5;1];

>> [RA,RB,n,X]=liezhu(A,b),H=det(A)

输出:

请注意:因为RA=RB,所以方程组有唯一解

RA =

4

RB =

4

n =

4

X =

0.0000

-1.0000

1.0000

1.0000

H =

-762.0001

第二题:

建立列主元高斯消去法m文件(题一中已有)

(1)输入:

>> format compact

>> A=[3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34];

>> b=[1;1;1];

>> [RA,RB,n,X]=liezhu(A,b),h=det(A),C=cond(A)

输出:

请注意:因为RA=RB,所以方程组有唯一解

RA =

3

RB =

n =

3

X =

1.0e+03 *

1.5926

-0.6319

-0.4936

h =

-0.0305

C =

3.0697e+04

(2)输入:

>> A=[3.00 6.03 1.99;1.27 4.16 -1.23;0.990 -4.81 9.34]; >> b=[1;1;1];

>> [RA,RB,n,X]=liezhu(A,b),h=det(A)

输出:

请注意:因为RA=RB,所以方程组有唯一解

RA =

3

RB =

3

n =

3

X =

119.5273

-47.1426

-36.8403

h =

-0.4070

第三题:

输入:

>> clear

>> A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];

>> b=[32 23 33 31]’;

>> dA=det(A),lamda=eig(A),Ac2=cond(A,2)

输出:

dA =

1.0000

lamda =

0.0102

0.8431

3.8581

30.2887

2.9841e+03

下面分析误差性态:

建立m文件:

function Acp=pjwc(A,jA,b,jb,p)

%Acp矩阵A的p条件数cond

%pjwc:p范数解的误差性态分析

%jA是A的近似矩阵jA=A+δA,jb=b+δb

Acp=cond(A,p);dA=det(A);X=A\b;

deltaA=jA-A;

pndA=norm(deltaA,p);deltab=jb-b;

pndb=norm(deltab,p);

if pndb>0

jX=A\jb;Pnb=norm(b,p);pnjx=norm(jX,p);deltaX=jX-X;

pnjdX=norm(deltaX,p);jxX=pnjdX/pnjX;

pnX=norm(X,p);xX=pnjdX/pnX;

pndb=norm(deltab,p);xAb=pndb/pnb;pnbj=norm(jb,p);xAbj=pndb/pnbj;

Xgxx=Acp*xAb;

end

if pndA>0

jX=jA\b;deltaX=jX-X;pnX=norm(X,p);

pnjdX=norm(deltaX,p);

pnjX=norm(jX,p);jxX=pnjdX/pnjX;xX=pnjdX/pnX;

pnjA=norm(jA,p);pnA=norm(A,p);

pndA=norm(deltaA,p);xAbj=pndA/pnjA;xAb=pndA/pnA;

Xgxx=Acp*xAb;

end

if (Acp>50)&(dA<0.1)

disp('请注意:AX=b是病态的,A的p条件数Acp,A的行列式值dA,解X,近似解jX,解的相对误差xX,解的相对误差估计Xgxx,b或A的相对误差xAb依次如下:')

Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'

else

disp('请注意:AX=b是良态的,A的p条件数Acp,A的行列式值dA,解X,近似解jX,解的相对误差xX,解的相对误差估计Xgxx,b或A的相对误差xAb依次如下:')

Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'

end

输入:

>> jA=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98];

>> jb=b;p=2;

>> Acp=pjwc(A,jA,b,jb,p)

输出:

请注意:AX=b是良态的,A的p条件数Acp,A的行列式值dA,解X,近似解jX,解的相对误差xX,解的相对误差估计Xgxx,b或A的相对误差xAb依次如下:

Acp =

2.9841e+03

dA =

1.0000

ans =

1.0000 1.0000 1.0000 1.0000 ans =

-9.5863 18.3741 -3.2258 3.5240 xX =

10.4661

jxX =

0.9842

Xgxx =

22.7396

xAb =

0.0076

xAbj =

0.0076

Acp =

2.9841e+03

第四题:

(1)输入:

建立m文件:

for n=2:6

a=hilb(n);

pnH(n-1)=cond(a,inf);

end

pnH

n=2:6;

plot(n,pnH);

数值分析(第五版)计算实习题第五章作业

可见条件数随着n的增大而急剧增大

(2)输入:

>> n=2;H=hilb(n);

>> x=(linspace(1,1,n))';

>> b=H*x;

>> [RA,RB,n,X]=gauss(H,b)

输出:

请注意:因为RA=RB,所以方程组有唯一解RA =

2

RB =

2

n =

2

X =

1.0000

1.0000

输入:>> r=b-H*X,deltax=X-x

输出:

r =

deltax =

1.0e-15 *

0.4441

-0.6661

输入:

>> n=3;H=hilb(n);

>> x=(linspace(1,1,n))';

>> b=H*x;

>> [RA,RB,n,X]=gauss(H,b)

>> r=b-H*X,deltax=X-x

输出:

X =

1.0000

1.0000

1.0000

r =

1.0e-15 *

0.2220

deltax =

1.0e-13 *

-0.0200

0.1221

-0.1255

>> n=4;H=hilb(n);

>> x=(linspace(1,1,n))'; >> b=H*x;

>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-x

X =

1.0000

1.0000

1.0000

1.0000

r =

1.0e-15 *

-0.4441

-0.1110

deltax =

1.0e-12 *

-0.0222

0.2485

-0.5980

0.3886

>> n=5;H=hilb(n);

>> x=(linspace(1,1,n))'; >> b=H*x;

>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-x

X =

1.0000

1.0000

1.0000

1.0000

1.0000

r =

1.0e-15 *

0.2220

0.1110

deltax =

1.0e-11 *

-0.0035

0.0524

-0.1937

0.2591

-0.1148

>> n=6;H=hilb(n);

>> x=(linspace(1,1,n))'; >> b=H*x;

>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-x

X =

1.0000

1.0000

1.0000

1.0000

1.0000

r =

1.0e-15 *

0.2220

0.1110

deltax =

1.0e-11 *

-0.0035

0.0524

-0.1937

0.2591

-0.1148

>> n=7;H=hilb(n);

>> x=(linspace(1,1,n))'; >> b=H*x;

>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-x

X =

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

r =

1.0e-15 *

0.2220

0.1110

deltax =

1.0e-09 *

-0.0008

0.0219

-0.1482

0.3854

-0.4254

0.1677

>> n=8;H=hilb(n);

>> x=(linspace(1,1,n))'; >> b=H*x;

>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-x

X =

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

r =

1.0e-15 *

-0.2220

-0.1110

-0.1110

deltax =

1.0e-06 *

-0.0000

0.0018

-0.0236

0.1279

-0.3442

0.4870

-0.3466

0.0978

>> n=9;H=hilb(n);

>> x=(linspace(1,1,n))'; >> b=H*x;

>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-x

X =

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

r =

1.0e-15 *

0.4441

0.2220

-0.2220

0.2220

0.2220

-0.1110

deltax =

1.0e-04 *

-0.0000

0.0002

-0.0028

0.0197

-0.0722

0.1471

-0.1687

0.1017

-0.0251

>> n=10;H=hilb(n);

>> x=(linspace(1,1,n))'; >> b=H*x;

>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-x

X =

1.0000

1.0000

1.0000

1.0000

0.9999

1.0003

0.9996

1.0004

0.9998

1.0000

r =

1.0e-15 *

0.4441

-0.2220

0.2220

-0.1110

0.1110

0.1110

deltax =

1.0e-03 *

-0.0000

0.0001

-0.0023

0.0205

-0.0974

0.2669

-0.4369

0.4214

-0.2209

0.0485