数值分析第五章答案.docx
《数值分析第五章答案.docx》由会员分享，可在线阅读，更多相关《数值分析第五章答案.docx（13页珍藏版）》请在冰豆网上搜索。
数值分析第五章答案
数值分析第五章答案
【篇一：
数值分析第五版计算实习题】
第二章
2-1
程序：
clear;clc;
x1=[0.20.40.60.81.0];
y1=[0.980.920.810.640.38];
n=length（y1）;
c=y1（:
）;
orj=2:
n%求差商
fori=n:
-1:
j
c（i）=（c（i）-c（i-1））/（x1（i）-x1（i-j+1））;
end
end
symsxdfd;
df
（1）=1;d
（1）=y1
（1）;
fori=2:
n%求牛顿差值多项式
df（i）=df（i-1）*（x-x1（i-1））;
d（i）=c（i）*df（i）;
end
disp（4次牛顿插值多项式）;
p4=vpa（collect（（sum（d）））,5）%p4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape（x1,y1,variational）;%调用三次样条函数
q=pp.coefs;
disp（三次样条函数）;
fori=1:
4
s=q（i,:
）*[（x-x1（i））^3;（x-x1（i））^2;（x-x1（i））;1];
s=vpa（collect（s）,5）
end
x2=0.2:
0.08:
1.08;
dot=[121112];
figure
ezplot（p4,[0.2,1.08]）;
holdon
y2=fnval（pp,x2）;
x=x2（dot）;
y3=eval（p4）;
y4=fnval（pp,x2（dot））;
plot（x2,y2,r,x2（dot）,y3,b*,x2（dot）,y4,co）;
title（4次牛顿插值及三次样条）;
结果如下：
4次牛顿插值多项式
p4=-0.52083*x^4+0.83333*x^3-1.1042*x^2+0.19167*x+0.98三次样条函数
x∈[0.2,0.4]时，s=-1.3393*x^3+0.80357*x^2-0.40714*x+1.04x∈[0.4,0.6]时，s=0.44643*x^3-1.3393*x^2+0.45*x+0.92571x∈[0.6,0.8]时，s=-1.6964*x^3+2.5179*x^2-1.8643*x+1.3886x∈[0.8,1.0]时，s=2.5893*x^3-7.7679*x^2+6.3643*x-0.80571输出图如下
2-3
（1）
程序：
clear;
clc;
x1=[01491625364964];
y1=[012345678];%插值点
n=length（y1）;
a=ones（n,2）;
a（:
2）=-x1;
c=1;
fori=1:
n
c=conv（c,a（i,:
））;
end
q=zeros（n,n）;
r=zeros（n,n+1）;
fori=1:
n
[q（i,:
）,r（i,:
）]=deconv（c,a（i,:
））;%wn+1/（x-xk）
end
dw=zeros（1,n）;
fori=1:
n
dw（i）=y1（i）/polyval（q（i,:
）,x1（i））;%系数
end
p=dw*q;
symsxl8;
fori=1:
n
l8（i）=p（n-i+1）*x^（i-1）;
end
disp（8次拉格朗日插值）;
l8=vpa（collect（（sum（l8）））,5）
xi=0:
64;
yi=polyval（p,xi）;
figure
plot（xi,yi,x1,y1,r*）;
holdon
title（8次拉格朗日插值）;
结果如下：
8次拉格朗日插值
l8=-3.2806e-10*x^8+6.7127e-8*x^7-5.4292e-6*x^6+0.00022297*x^5-0.0049807*x^4+0.060429*x^3-0.38141*x^2+1.3257*x
输出图如下：
第五章
4-1（3）
程序：
clc;
clear;
y=@（x）sqrt（x）.*log（x）;
a=0;b=1;tol=1e-4;
p=quad（y,a,b,tol）;
fprintf（采用自适应辛普森积分结果为：
%d\n,p）;
结果如下：
采用自适应辛普森积分结果为：
-4.439756e-01
第九章
9-1
（a）程序：
clc;
clear;
a=1;b=2;%定义域
h=0.05;%步长
n=（b-a）/h;
y0=1;%初值
f=@（x,y）1/x^2-y/x;%微分函数
xn=linspace（a,b,n+1）;%将定义域分为n等份yn=zeros（1,n）;%结果矩阵
yn
（1）=y0;%赋初值
%以下根据改进欧拉公式求解
fori=1:
n
xn=xn（i）;
xnn=xn（i+1）;
yn=yn（i）;
yp=yn+h*f（xn,yn）;
yc=yn+h*f（xnn,yp）;
yn=（yp+yc）/2;
yn（i+1）=yn;
end
xn=yn;
%以下根据经典四阶r-k法公式求解
fori=1:
n
xn=xn（i）;
yn=yn（i）;
k1=f（xn,yn）;
k2=f（xn+h/2,yn+h/2*k1）;
k3=f（xn+h/2,yn+h/2*k2）;
k4=f（xn+h,yn+h*k3）;
yn=yn+h/6*（k1+2*k2+2*k3+k4）;
yn（i+1）=yn;
end
disp（改进欧拉法四阶经典r-k法）;disp（[xnyn]）
结果如下：
改进欧拉法四阶经典r-k法11
0.998870.99885
0.995770.9978
0.991140.99694
0.985320.99634
0.978570.99603
0.971110.99606
0.963110.99645
0.95470.99723
0.945980.99841
0.937051
0.927981.002
0.918831.0044
0.909641.0073
0.900451.0106
0.891291.0143
0.882181.0184
0.873151.0229
0.864211.0278
0.855381.0331
0.846651.0388
（b）程序：
clc;
clear;
a=0;b=1;%定义域
h=[0.10.0250.01];%步长
y0=1/3;%初值
f=@（x,y）-50*y+50*x^2+2*x;%微分函数xi=linspace（a,b,11）;
y=1/3*exp（-50*xi）+xi.^2;%准确解ym=zeros（1,11）;
forj=1:
3
【篇二：
数值分析（第五版）计算实习题第五章作业】
题：
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;
ifij
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;
ifzhicha0
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=rbn，所以方程组有无穷多解）
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=
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
2.9841e+03
下面分析误差性态：
建立m文件：
functionacp=pjwc（a,ja,b,jb,p）
%acp矩阵a的p条件数cond
%pjwc：
p范数解的误差性态分析
acp=cond（a,p）;da=det（a）;x=a\b;
deltaa=ja-a;
pnda=norm（deltaa,p）;deltab=jb-b;
pndb=norm（deltab,p）;
ifpndb0
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
ifpnda0
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（acp50）（da0.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=
【篇三：
李庆扬-数值分析第五版第5章与第7章习题答案】
升级会员 