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

数值分析(第五版)计算实习题 matlab编程

数值分析第五章

第一题：

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

数值分析(第五版)计算实习题 matlab编程

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;

数值分析(第五版)计算实习题 matlab编程

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<n，所以方程组有无穷多解')

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

>> [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 =

数值分析(第五版)计算实习题 matlab编程

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

数值分析(第五版)计算实习题 matlab编程

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 =…… 此处隐藏：4582字，全部文档内容请下载后查看。喜欢就下载吧 ……