数值分析第五版计算实习题第五章作业.docx
《数值分析第五版计算实习题第五章作业.docx》由会员分享，可在线阅读，更多相关《数值分析第五版计算实习题第五章作业.docx（18页珍藏版）》请在冰豆网上搜索。
数值分析第五版计算实习题第五章作业
数值分析第五章
第一题：
LU分解法：
建立m文件
functionh1=zhijieLU（A,b）%h1各阶主子式的行列式值
[nn]=size（A）;RA=rank（A）;
ifRA~=n
disp（'请注意：
因为A的n阶行列式h1等于零，所以A不能进行LU分解。
A的秩RA如下：
'）
RA,h1=det（A）;
return
end
ifRA==n
forp=1:
n
h（p）=det（A（1:
p,1:
p））;
end
h1=h（1:
n）;
fori=1:
n
ifh（1,i）==0
disp（'请注意：
因为A的r阶主子式等于零，所以A不能进行LU分解。
A的秩RA和各阶顺序主子式h1依次如下：
'）
h1;RA
return
end
end
ifh（1,i）~=0
disp（'请注意：
因为A的r阶主子式都不等于零，所以A能进行LU分解。
A的秩RA和各阶顺序主子式h1依次如下：
'）
forj=1:
n
U（1,j）=A（1,j）;
end
fork=2:
n
fori=2:
n
forj=2:
n
L（1,1）=1;L（i,i）=1;
ifi>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-701;-32.09999962;5-15-1;2102];
>>b=[8;5.900001;5;1];
>>h1=zhijieLU（A,b）
输出：
请注意：
因为A的r阶主子式都不等于零，所以A能进行LU分解。
A的秩RA和各阶顺序主子式h1依次如下：
RA=
4
U=
10.0000-7.000001.0000
02.10006.00002.3000
00-2.1429-4.2381
0-0.0000012.7333
L=
1.0000000
-0.30001.000000
0.50001.19051.0000-0.0000
0.20001.14293.20001.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=[Ab];n=length（b）;RA=rank（A）;RB=rank（B）;zhicha=RB-RA;
ifzhicha>0
disp（'请注意：
因为RA~=RB，所以方程组无解'）
return
warningoffMATLAB:
return_outside_of_loop
end
ifRA==RB
ifRA==n
disp（'请注意：
因为RA=RB，所以方程组有唯一解'）
X=zeros（n,1）;C=zeros（1,n+1）;
forp=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;
fork=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）;
forq=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=RBend
end
输入：
>>A=[10-701;-32.09999962;5-15-1;2102];
>>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）输入：
>>formatcompact
>>A=[3.016.031.99;1.274.16-1.23;0.987-4.819.34];
>>b=[1;1;1];
>>[RA,RB,n,X]=liezhu（A,b）,h=det（A）,C=cond（A）
输出：
请注意：
因为RA=RB，所以方程组有唯一解
RA=
3
RB=
3
n=
3
X=
1.0e+03*
1.5926
-0.6319
-0.4936
h=
-0.0305
C=
3.0697e+04
（2）输入：
>>A=[3.006.031.99;1.274.16-1.23;0.990-4.819.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=[10787;7565;86109;75910];
>>b=[32233331]’;
>>dA=det（A）,lamda=eig（A）,Ac2=cond（A,2）
输出：
dA=
1.0000
lamda=
0.0102
0.8431
3.8581
30.2887
Ac2=
2.9841e+03
下面分析误差性态：
建立m文件：
functionAcp=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）;
ifpndb>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
ifpndA>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=[1078.17.2;7.085.0465;85.989.899;6.99599.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.00001.00001.00001.0000
ans=
-9.586318.3741-3.22583.5240
xX=
10.4661
jxX=
0.9842
Xgxx=
22.7396
xAb=
0.0076
xAbj=
0.0076
Acp=
2.9841e+03
第四题：
（1）输入：
建立m文件:
forn=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=
0
0
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
0
0
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
-0.1110
0
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
0.2220
0
0
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
0.2220
0
0
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
0
0.2220
0
0
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
0
-0.2220
0
0
-0.1110
-0.1110
0
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
0.2220
0.2220
0
-0.1110
0
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
0
-0.2220
0.2220
0
0
-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
升级会员 