数值分析第五版计算实习题第五章作业教学资料.docx
- 文档编号:4742541
- 上传时间:2022-12-08
- 格式:DOCX
- 页数:18
- 大小:24.67KB
数值分析第五版计算实习题第五章作业教学资料.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=RB end 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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 第五 计算 实习 作业 教学 资料