矩阵与数值分析作业.docx
《矩阵与数值分析作业.docx》由会员分享，可在线阅读，更多相关《矩阵与数值分析作业.docx（12页珍藏版）》请在冰豆网上搜索。
矩阵与数值分析作业
2011级工科硕士研究生
《矩阵与数值分析》课程数值实验题目
一、对于数列
，有如下两种生成方式
1、首项为
，递推公式为
；
2、前两项为
，递推公式为
；
给出利用上述两种递推公式生成的序列的第50项。
解：
matlab编程
1、
c1
（1）=1;
fort1=2:
50
c1（t1）=c1（t1-1）/3;
end
c1（50）
运算结果：
ans=4.1788e-024
2、
c2
（1）=1;c2
（2）=1/3;
fort2=3:
50
c2（t2）=c2（t2-1）*10/3-c2（t2-2）;
end
c2（50）
运算结果：
ans=-4.9661e+006
二、利用迭代格式
及Aitken加速后的新迭代格式求方程
在
内的根
解：
matlab编程
k=1;x
（1）=1.3;f=1;
while（abs（f）>1.0e-6）
x（k+1）=sqrt（10/（x（k）+4））;
k=k+1;
f=x（k）^3+4*x（k）^2-10;
end
x,f
clear
k=1;x
（1）=1.3;f=1;
while（abs（f）>1.0e-6）
y（k）=sqrt（10/（x（k）+4））;
z（k）=sqrt（10/（y（k）+4））;
x（k+1）=z（k）-（z（k）-y（k））^2/（z（k）-2*y（k）+x（k））;
k=k+1;
f=x（k）^3+4*x（k）^2-10;
end
x,f
运算结果：
x=1.30001.37351.36411.36531.36521.36511.36511.3651
f=5.8603e-007
x=1.30001.36511.3651
f=2.5046e-012
分析：
原迭代格式运算8次，达到精度要求；而是用Aitken加速后，只用3次迭代，即Aitken加速可以明显加快迭代速度。
三、解线性方程组
1．分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组
迭代法计算停止的条件为：
．
解：
matlab编程
1
Xj=[0;0;0;0];Xg=[0;0;0;0];
A=[621-2;
250-2;
-2085;
1327];
b=[4;7;-1;0];
D=diag（diag（A））;
L=-tril（A）+D;
U=-triu（A）+D;
Bj=inv（D）*（L+U）;
f=inv（D）*b;
k=1;flag=1;
Xj（:
k+1）=Bj*Xj（:
k）+f;
whileflag==1
ifabs（Xj（:
k+1）-Xj（:
k））<1e-6
flag=0;
else
k=k+1;
Xj（:
k+1）=Bj*Xj（:
k）+f;
end
end
Xj
Bg=inv（D-L）*U;
fg=inv（D-L）*b;
k=1;flag=1;
Xg（:
k+1）=Bg*Xg（:
k）+fg;
whileflag==1
ifabs（Xg（:
k+1）-Xg（:
k））<1e-6
flag=0;
else
k=k+1;
Xg（:
k+1）=Bg*Xg（:
k）+fg;
end
end
Xg
运算结果：
Xj=Columns1through11
00.66680.22080.06220.08380.05670.05370.05350.05220.05220.0521
01.40001.13331.04781.16351.14421.14781.15151.15051.15091.1509
0-0.12500.04160.34240.22120.24330.24780.24320.24480.24470.2446
00-0.6593-0.5292-0.5557-0.5738-0.5681-0.5704-0.5706-0.5705-0.5706
Columns12through16
0.05210.05210.05210.05200.0520
1.15091.15091.15091.15091.1509
0.24450.24450.24450.24450.2445
-0.5707-0.5707-0.5707-0.5707-0.5707
Xg=
00.66670.08430.05510.05230.05210.05210.05200.0520
01.13321.12911.14901.15081.15081.15091.15091.1509
00.04170.26660.24640.24480.24460.24460.24460.2446
0-0.5928-0.5721-0.5707-0.5706-0.5706-0.5706-0.5706-0.5706
2.用Gauss列主元消去法、QR方法求解如下方程组：
2、
A=[2212;
413-1;
-4-201;
2323];
b=[1;2;1;0];
Ab=[Ab];
forj=1:
3
k=j;
max（j）=Ab（j,j）;
fori=（j+1）:
4
ifmax（j）k=i;
max（j）=Ab（i,j）;
end
end
ifk~=j
t=Ab（j,:
）;
Ab（j,:
）=Ab（k,:
）;
Ab（k,:
）=t;
end
fori=（j+1）:
4
Ab（i,:
）=Ab（i,:
）-Ab（j,:
）*Ab（i,j）/Ab（j,j）;
end
end
fori=4:
-1:
1
Abt=Ab（i,5）;
forj=4:
-1:
（i+1）
Abt=Abt-x（j）*Ab（i,j）;
end
x（i）=Abt/Ab（i,i）;
end
Xg=x'
[Q,R]=qr（A）;
Xqr=R\Q'*b
运算结果：
Xg=
1.5417
-2.7500
0.0832
1.6667
Xqr=
1.5417
-2.7500
0.0832
1.6667
分析：
G-S迭代法与Jacob迭代法相比结果相同，但减少很多迭代步数。
四、已知一组数据点
，编写一程序求解三次样条插值函数
满足
并针对下面一组具体实验数据
0.25
0.3
0.39
0.45
0.53
0.5000
0.5477
0.6245
0.6708
0.7280
求解，其中边界条件为
.
解：
n=5;
X=[0.250.30.390.450.53];
Y=[0.50000.54770.62450.67080.7280];
f00=0;
fnn=0;
fork=1:
（n-1）
h（k）=X（k+1）-X（k）;
end
g
（1）=3*（Y
（2）-Y
（1））/h
（1）-f00*h
（1）/2;
g（n）=3*（Y（n）-Y（n-1））/h（n-1）+fnn*h（n-1）/2;
fork=2:
（n-1）
v（k）=h（k）/（h（k）+h（k-1））;
u（k）=h（k-1）/（h（k）+h（k-1））;
g（k）=3*u（k）*（Y（k+1）-Y（k））/h（k）+3*v（k）*（Y（k）-Y（k-1））/h（k-1）;
end
u
（1）=1;
v
（1）=[];v（n-1）=1;
A=2*eye（n）+diag（u,1）+diag（v,-1）;
m=A\g';
symsxSx;
fork=1:
（n-1）
fprintf（'%gSx=（h（k）+2*（x-X（k）））*（x-X（k+1））^2*Y（k）/h（k）^3+...
（h（k）-2*（x-X（k+1）））*（x-X（k））^2*Y（k+1）/h（k）^3+...
（x-X（k））*（x-X（k+1））^2*m（k）/h（k）^2+...
（x-X（k+1））*（x-X（k））^2*m（k+1）/h（k）^2;
Sx=vpa（expand（Sx）,4）
End
运算结果：
0.25Sx=
4.699*x^2-0.2051*x+0.3555-6.265*x^3
0.3Sx=
-2.633*x^2+1.995*x+0.1355+1.881*x^3
0.39Sx=
0.1064*x^2+0.9261*x+0.2744-0.4600*x^3
0.45Sx=
-3.409*x^2+2.508*x+0.3710e-1+2.144*x^3
分析：
三次样条函数具有良好的光滑性和几何光顺性，同时收敛性和稳定性也较好。
五、编写程序构造区间
上的以等分结点为插值结点的Newton插值公式，假设结点数为
（包括两个端点），给定相应的函数值，插值区间和等分的份数，该程序能快速计算出相应的插值公式。
以
，
为例计算其对应的插值公式，分别取不同的
值并画出原函数的图像以及插值函数的图像，观察当
增大时的逼近效果.
解：
matlab编程
functionLagrange
clc;clear;closeall;
fori=1:
3;
ifi==1
N=5;
elseifi==2
N=8;
else
N=10;
end
f=inline（'1/（1+25*x^2）'）;
x1=zeros（1,N+1）;
a=-1;
b=1;
fori=1:
N+1
x1（i）=a+（i-1）*（b-a）/N;
y1（i）=feval（f,x1（i））;
end
symsx
ff=0;
fori=1:
N+1
f=1;
forj=1:
i-1
f=f*（x-x1（j））/（x1（i）-x1（j））;
end
forj=i+1:
N+1
f=f*（x-x1（j））/（x1（i）-x1（j））;
end
ff=f*y1（i）+ff;
f=1;
end
ff=collect（ff,x）;
ff=vpa（ff,4）;
y=ff;
p=ezplot（y,[a,b]）;grid
YLIM（[-0.10.6]）;
ifN==5
set（p,'Color','black'）;
set（p,'LineStyle','--'）;
lagrange_5=y
elseifN==8
set（p,'Color','r'）;
set（p,'LineStyle','--'）;
lagrange_8=y
else
set（p,'Color','b'）;
set（p,'LineStyle','--'）
lagrange_10=y
end
holdon;
xlabel（'x'）;ylabel（'y'）;
title（'y=p（x）'）;holdon;
end
Lag_Cheb（）;
x=-1:
0.01:
1;
y=1./（1+25*x.^2）;
acu=plot（x,y）;gridon;holdon
set（acu,'Color','m'）;
set（acu,'LineStyle','-'）;
legend（'N=5','N=8','N=10','Cheb,N=10','¾«È·Öµ'）;
functionLag_Cheb（）
f=inline（'1/（1+25*x^2）'）;
N=10;
x1=zeros（1,N+1）;
a=-1;
b=1;
fori=1:
N+1
x1（i）=cos（（2*i-1）*pi/（2*N））;
y1（i）=feval（f,x1（i））;
end
symsx
ff=0;
fori=1:
N+1
f=1;
forj=1:
i-1
f=f*（x-x1（j））/（x1（i）-x1（j））;
end
forj=i+1:
N+1
f=f*（x-x1（j））/（x1（i）-x1（j））;
end
ff=f*y1（i）+ff;
f=1;
end
ff=collect（ff,x）;
ff=vpa（ff,4）;
yy=ff;
ff=collect（ff,x）;
yy=ff;
lagrange_chebshev_10=yy
cheb=ezplot（yy,[a,b]）;gridon
YLIM（[-0.10.6]）;
set（cheb,'Color','g'）;
set（cheb,'LineStyle',':
'）;
结果输出
lagrange_5=0.5673+1.202*x^4-1.731*x^2
lagrange_8=1.+53.69*x^8-102.8*x^6+61.37*x^4-13.20*x^2
lagrange_10=1.-220.9*x^10+494.9*x^8-381.4*x^6+123.4*x^4-16.86*x^2
lagrange_chebshev_10=0.7413-5.359*x^10-0.4e-2*x^9+18.96*x^8-0.1321*x^7-
25.78*x^6+0.10*x^5+16.81*x^4+0.5e-2*x^3-5.336*x^2+0.5288e-3*x
分析：
（1）拉格朗日插值法随着插值节点增多而更加接近真实曲线，但是还是会在区间
［-1，-0.4］及［0.4，1］上，由于高次振荡产生Runge现象。
（2）分段插值法则优于前者，因为分段插值在每个小区间上都可以最大限度的满足边界条件，但光滑度较低。
升级会员 