矩阵与数值分析.docx
《矩阵与数值分析.docx》由会员分享，可在线阅读，更多相关《矩阵与数值分析.docx（25页珍藏版）》请在冰点文库上搜索。
矩阵与数值分析
2013级工科硕士研究生
《矩阵与数值分析》课程数值实验题目
一、设
，分别编制从小到大和从大到小的顺序程序分别计算
并指出两种方法计算结果的有效位数。
Matlab程序如下：
function[si,sd]=S（N）
formatlong;
si=0;sd=0;
forj=N:
-1:
2
si=1.0e6/（j^2-1）+si;
end
forj=2:
N
sd=1.0e6/（j^2-1）+sd;
end
end
在matlab命令窗口中输入：
[si,sd]=S（10000）
运行结果：
si=7.499000049995000e+005
sd=7.499000049994994e+005
在matlab命令窗口中输入：
[si,sd]=S（1000000）
运行结果：
si=7.499990000005000e+005
sd=7.499990000005200e+0051
结果分析：
si为从大到小的顺序求和的值，sd为从小到大的顺序求和的值。
当N分别为10000和1000000时，si分别为7.499000049995000e+005和7.499990000005000e+005，可以看出这两个数的有效值均为13位；而sd分别为7.499000049994994e+005和7.499990000005200e+005，这两个数的有效值均为16位。
这就出现了我们在矩阵理论课上所学的“大数吃小数”的问题。
为了使结果更为精确我们必须避免在四则运算中出现“大数吃小数”的情况，应该按从小到大的顺序进行求和。
二、解线性方程组
1．分别利用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组
，其中常向量为
维随机生成的列向量，系数矩阵
具有如下形式
，
其中
为
阶矩阵，
为
阶单位矩阵，迭代法计算停止的条件为：
，给出
时的不同迭代步数．
求解系数矩阵A和向量b的Matlab程序如下：
function[A,b,x0]=jz（n）
fori=1:
n-1%此处n赋值即可，如n=100
forj=1:
n-1
if（i==j）
T（i,j）=2;
end
if（i==（j+1））
T（i,j）=-1;
end
if（i==（j-1））
T（i,j）=-1;
end
end
end
I=eye（n-1）;
k=1;
formm=1:
（n-1）
A（k:
（k+n-2）,k:
（k+n-2））=T+2*I;
k=k+n-1;
end
k=1;
forxx=1:
（n-2）
A（k:
（k+n-2）,（k+n-1）:
（k+2*n-3））=-I;
A（（k+n-1）:
（k+2*n-3）,k:
（k+n-2））=-I;
k=k+n-1;
end
b=randn（（n-1）^2,1）;
x0=zeros（（n-1）^2,1）;
Jacobi迭代法Matlab程序如下：
functionn=jacobi（A,b,x0）
eps=1.0e-10;
M=10000;
D=diag（diag（A））;%求A的对角矩阵
L=-tril（A,-1）;%求A的下三角阵
U=-triu（A,1）;%求A的上三角阵
B=D\（L+U）;
f=D\b;
x=B*x0+f;
n=1;%迭代次数
whilenorm（x-x0）>=eps
x0=x;
x=B*x0+f;
n=n+1;
if（n>=M）
disp（'Warning:
迭代次数太多，可能不收敛！
'）;
return;
end
end
Gauss-Seidel迭代法Matlab程序如下：
functionn=gauseidel（A,b,x0）
eps=1.0e-10;
M=10000;
D=diag（diag（A））;%求A的对角矩阵
L=-tril（A,-1）;%求A的下三角阵
U=-triu（A,1）;%求A的上三角阵
G=（D-L）\U;
f=（D-L）\b;
x=G*x0+f;
n=1;%迭代次数
whilenorm（x-x0）>=eps
x0=x;
x=G*x0+f;
n=n+1;
if（n>=M）
disp（'Warning:
迭代次数太多，可能不收敛！
'）;
return;
end
end
在matlab命令窗口中输入：
[A,b,x0]=jz（10）;
J10=jacobi（A,b,x0）
G10=gauseidel（A,b,x0）
[A,b,x0]=jz（20）;
J20=jacobi（A,b,x0）
G20=gauseidel（A,b,x0）
[A,b,x0]=jz（30）;
J30=jacobi（A,b,x0）
G30=gauseidel（A,b,x0）
运行结果：
J10=436；G10=214
J20=1810；G20=913
J30=3990；G30=2001
结果分析：
N=10且M=1000时，Jacobi迭代法和Gauss—seidel迭代法的迭代次数分别为436和214；N=20且M=1000时，Jacobi迭代法和Gauss—seidel迭代法的迭代次数分别为1810和913；N=30且M=1000时，Jacobi和Gauss-seidel算法的迭代次数分别为3990和2001次。
从以上结果可知在该问题下Jacobi算法的收敛性没有Gauss-seidel算法的收敛性好，原因在于Jacobi算法迭代过程中，对已算出的信息未加充分利用，一般来说，后面计算的值要比前面的计算值要精确些，而Gauss-seidel算法则充分利用了已求出来的信息，所以此算法的收敛性更为好一些。
2.用Gauss列主元消去法、QR方法求解如下方程组：
Gauss列主元消去法Matlab程序如下：
function[x,XA]=GaussXQLineMain（A,b）
N=size（A）;
n=N
（1）;
index=0;
fori=1:
（n-1）
me=max（abs（A（1:
n,i）））;%选取列主元
fork=i:
n
if（abs（A（k,i））==me）
index=k;%保存列主元所在的行
break;
end
end
temp=A（i,1:
n）;
A（i,1:
n）=A（index,1:
n）;
A（index,1:
n）=temp;
bb=b（index）;
b（index）=b（i）;
b（i）=bb;%交换主行
forj=（i+1）:
n
if（A（i,i）==0）
disp（'对角元素为0！
'）;
return;
end
l=A（j,i）;
m=A（i,i）;
A（j,1:
n）=A（j,1:
n）-l*A（i,1:
n）/m;
b（j）=b（j）-l*b（i）/m;%消元
end
end
fori=n:
-1:
1
if（is=A（i,（i+1）:
n）*x（（i+1）:
n,1）;
else
s=0;
end
x（i,1）=（b（i）-s）/A（i,i）;
end
XA=A;
QR法Matlab程序如下：
functionx=qrxq（A,b）
N=size（A）;
n=N
（1）;
B=A;%保存系数矩阵
A（1:
n,1）=A（1:
n,1）/norm（A（1:
n,1））;%将A的第一列正规化
fori=2:
n
forj=1:
（i-1）
A（1:
n,i）=A（1:
n,i）-dot（A（1:
n,i）,A（1:
n,j））*A（1:
n,j）;
%使A的第i列与前面所有的列正交
end
A（1:
n,i）=A（1:
n,i）/norm（A（1:
n,i））;
%将A的第i列正规化
end
Q=A;%分解出来的正交矩阵
R=transpose（Q）*B;
bb=transpose（Q）*b;
fori=n:
-1:
1
if（is=R（i,（i+1）:
n）*x（（i+1）:
n,1）;
else
s=0;
end
x（i,1）=（bb（i）-s）/R（i,i）;
end
在matlab命令窗口中输入：
A=[12-11;2503;1792;8-1-21];
b=[0;4;12;-8];
x_G=GaussXQLineMain（A,b）
x_QR=qrxq（A,b）
运行结果：
x_G=-1.00000000000000
-0.00000000000000
1.00000000000000
2.00000000000000
x_QR=-1.00000000000000
-0.00000000000001
1.00000000000000
2.00000000000002
三、非线性方程的迭代解法
1．求方程
的根，迭代停止的条件为：
；
使用弦截法求解此方程的根，弦截法的Matlab程序如下：
functionroot=Secant（f,a,b）
%f:
方程;
%a:
区间左端点;
%b:
区间右端点;
%root:
求得的根;
eps=1.0e-10;
f1=subs（sym（f）,findsym（sym（f））,a）;
f2=subs（sym（f）,findsym（sym（f））,b）;
if（f1==0）
root=a;
end
if（f2==0）
root=b;
end
if（f1*f2>0）
disp（'两端点函数值乘积大于0!
'）;
return;
else
tol=1;
fa=subs（sym（f）,findsym（sym（f））,a）;
fb=subs（sym（f）,findsym（sym（f））,b）;
root=a-（b-a）*fa/（fb-fa）;
while（tol>eps）
r1=root;
fx=subs（sym（f）,findsym（sym（f））,r1）;
s=fx*fa;
if（s==0）
root=r1;
else
if（s>0）
root=b-（r1-b）*fb/（fx-fb）;
else
root=a-（r1-a）*fa/（fx-fa）;
end
end
tol=abs（root-r1）;
end
end
利用matlab的绘图功能可以确定求根区间为（1,3）
在matlab命令窗口中输入：
root=Secant（'exp（x）+2*x^2+2*sin（x）-log（x）-16',1,3）
运行结果：
root=1.96292268316446
2．利用Newton迭代法求多项式
的所有实零点，注意重根的问题。
Newton迭代法的Matlab程序如下：
x=-3;d=-1;
whileabs（d-x）>0.5*10^（-8）%在区间（-3，-1）
d=x;
x=x-（x^4-3*x^3-3*x^2+11*x-6）/（4*x^3-9*x^2-6*x+11）;%牛顿法
end
x0=d
x=0;d=2;
whileabs（d-x）>0.5*10^（-8）%在区间（0，2）的重根
d=x;
x=x-2*（x^4-3*x^3-3*x^2+11*x-6）/（4*x^3-9*x^2-6*x+11）;
end
x1=d
x=2.5;d=4;
whileabs（d-x）>0.5*10^（-8）%在区间（2.5，4）
d=x;
x=x-（x^4-3*x^3-3*x^2+11*x-6）/（4*x^3-9*x^2-6*x+11）;
end
x2=d
运行结果：
x0=-2.00000000151233
x1=1.00000000065769
x2=3.00000000000002
结果分析：
由以上结果可知，newton迭代法所得结果与初始值的选取密切相关。
不同的初始值会产生不同的结果，甚至会影响牛顿迭代法收敛性。
四、数值积分
用三点Gauss型求积公式计算积分
n点Gauss型求积公式matlab程序如下：
functionI=question4（f,a,b,n）
%f为求积方程;
%a,b分别为积分上下限;
%n为Gauss点个数;
[XK,AK]=grule（n）;
XK1=（（b-a）/2）*XK+（（a+b）/2）;
I=（（b-a）/2）*sum（AK.*subs（f,findsym（f）,XK1））;
function[bp,wf]=grule（n）
%[bp,wf]=grule（n）
bp=zeros（n,1）;wf=bp;iter=2;m=fix（（n+1）/2）;e1=n*（n+1）;
mm=4*m-1;t=（pi/（4*n+2））*（3:
4:
mm）;nn=（1-（1-1/n）/（8*n*n））;
xo=nn*cos（t）;
forj=1:
iter
pkm1=1;pk=xo;
fork=2:
n
t1=xo.*pk;pkp1=t1-pkm1-（t1-pkm1）/k+t1;
pkm1=pk;pk=pkp1;
end
den=1.-xo.*xo;d1=n*（pkm1-xo.*pk）;dpn=d1./den;
d2pn=（2.*xo.*dpn-e1.*pk）./den;
d3pn=（4*xo.*d2pn+（2-e1）.*dpn）./den;
d4pn=（6*xo.*d3pn+（6-e1）.*d2pn）./den;
u=pk./dpn;v=d2pn./dpn;
h=-u.*（1+（.5*u）.*（v+u.*（v.*v-u.*d3pn./（3*dpn））））;
p=pk+h.*（dpn+（.5*h）.*（d2pn+（h/3）.*（d3pn+.25*h.*d4pn）））;
dp=dpn+h.*（d2pn+（.5*h）.*（d3pn+h.*d4pn/3））;
h=h-p./dp;xo=xo+h;
end
bp=-xo-h;
fx=d1-h.*e1.*（pk+（h/2）.*（dpn+（h/3）.*（...
d2pn+（h/4）.*（d3pn+（.2*h）.*d4pn））））;
wf=2*（1-bp.^2）./（fx.*fx）;
if（m+m）>n,bp（m）=0;end
if~（（m+m）==n）,m=m-1;end
jj=1:
m;n1j=（n+1-jj）;bp（n1j）=-bp（jj）;wf（n1j）=wf（jj）;
%end
在matlab命令窗口中输入：
symsx;
f=exp（-x）;
a=0;b=1;n=3;
I=question4（f,a,b,n）
运行结果：
I=0.632120255664068
五、插值与逼近
1．给定
上的函数
，请做如下的插值逼近：
（1）构造等距节点分别取
，
，
的Lagrange插值多项式；
（2）取Chebyshev多项式
的零点：
，
作插值节点构造
的插值多项式
和上述的插值多项式均要求画出曲线图形（用不同的线型或颜色表示不同的曲线）。
（1）Lagrange插值matlab程序如下：
functionf=L1（n）
symsts;
q=1/（1+s^2）;
fori=1:
n
x（i）=-1+2*（i-1）/（n-1）;
y（i）=subs（q,s,x（i））;
end
f=0.0;
for（i=1:
n）
l=y（i）;
for（j=1:
i-1）
l=l*（t-x（j））/（x（i）-x（j））;
end;
for（j=i+1:
n）
l=l*（t-x（j））/（x（i）-x（j））;%计算拉格朗日基函数
end;
f=f+l;%计算拉格朗日插值函数
f=collect（f）;%化简
end
X=-1:
0.1:
1;
Y=subs（q,s,X）;
F=subs（f,t,X）;
plot（X,Y,'g',X,F,'r'）;
f=collect（f）;%展开
digits（5）;
f=vpa（f）;
1、在matlab命令窗口中输入：
f=L1（5）
运行结果
2、在matlab命令窗口中输入：
f=L1（8）
运行结果
3、在matlab命令窗口中输入：
f=L1（10）
运行结果：
结果分析：
由图像可以看出n=5时拟合曲线与原始曲线的误差比较大，而n=8,10时拟合曲线与原始曲线几乎重合在一起。
可见节点个越多，拉格朗日插值法所得的多项式曲线与原始曲线越接近，误差就越小。
（2）取Chebyshev多项式
的零点为插值点,此时Lagrange插值matlab程序主体基本不变，变的只是插值点，其程序如下：
functionf=L2（n）
symsts;
q=1/（1+s^2）;
fori=1:
n
x（i）=cos（（2*i-1）*pi/2/n）;
y（i）=subs（q,s,x（i））;
end
f=0.0;
for（i=1:
n）
l=y（i）;
for（j=1:
i-1）
l=l*（t-x（j））/（x（i）-x（j））;
end;
for（j=i+1:
n）
l=l*（t-x（j））/（x（i）-x（j））;%计算拉格朗日基函数
end;
f=f+l;%计算拉格朗日插值函数
end
X=-1:
0.1:
1;%绘图
Y=subs（q,s,X）;
F=subs（f,t,X）;
plot（X,Y,'k',X,F,'rp'）;
f=collect（f）;%展开
digits（5）;
f=vpa（f）;
在matlab命令窗口中输入：
f=L2（10）
运行结果：
结果分析：
有图像比较可以看出，N=10时，chebyshev零为拉格朗日插值时，其拟合曲线与原始曲线几乎重合，精度较高。
2．已知函数值
3
4
6
6
0
2
和边界条件：
． 求三次样条插值函数
并画出其图形．
三次样条插值matlab程序如下：
functionThrSample1（x,y,y_1,y_N）
symst;
f=0.0;
if（length（x）==length（y））
n=length（x）;
else
disp（'x和y的维数不相等！
'）;
return;
end%维数检查
A=diag（2*ones（1,n））;%求解m的系数矩阵
u=zeros（n-2,1）;
lamda=zeros（n-1,1）;
c=zeros（n,1）;
fori=2:
n-1
u（i-1）=（x（i）-x（i-1））/（x（i+1）-x（i-1））;
lamda（i）=（x（i+1）-x（i））/（x（i+1）-x（i-1））;
c（i）=3*lamda（i）*（y（i）-y（i-1））/（x（i）-x（i-1））+...
3*u（i-1）*（y（i+1）-y（i））/（x（i+1）-x（i））;
A（i,i+1）=u（i-1）;
A（i,i-1）=lamda（i）;%形成系数矩阵及向量c
end
c
（1）=2*y_1;
c（n）=2*y_N;
m=followup（A,c）;%用追赶法求解方程组
j=1;
whilej<=n-1
h=x（j+1）-x（j）;
f=y（j）*（2*（t-x（j））+h）*（t-x（j+1））^2/h^3+...
y（j+1）*（2*（x（j+1）-t）+h）*（t-x（j））^2/h^3+...
m（j）*（t-x（j））*（x（j+1）-t）^2/h^2-...
m（j+1）*（x（j+1）-t）*（t-x（j））^2/h^2;
f=collect（f）;
fprintf（'当x属于[%d,%d],f=\n',x（j）,x（j+1））;
disp（f）;
X=x（j）:
0.1:
x（j+1）;%绘图
F=subs（f,t,X）;
plot（X,F,'g'）;
holdon;
j=j+1;
end
functionx=followup（A,b）
n=rank（A）;
for（i=1:
n）
if（A（i,i）==0）
disp（'Error:
对角有元素为0！
'）;
return;
end
end;
d=ones（n,1）;
a=ones（n-1,1）;
c=ones（n-1）;
for（i=1:
n-1）
a（i,1）=A（i+1,i）;
c（i,1）=A（i,i+1）;
d（i,1）=A（i,i）;
end
d（n,1）=A（n,n）;
for（i=2:
n）
d（i,1）=d（i,1）-（a（i-1,1）/d（i-1,1））*c（i-1,1）;
b（i,1）=b（i,1）-（a（i-1,1）/d（i-1,1））*b（i-1,1）;
end
x（n,1）=b（n,1）/d（n,1）;
for（i=（n-1）:
-1:
1）
x（i,1）=（b（i,1）-c（i,1）*x（i+1,1））/d（i,1）;
end
在matlab命令窗口中输入：
x=[346];
y=[602];
y_1=1;
y_N=-1;
ThrSample1（x,y,y_1,y_N）
运行结果：
结果分析：
由图像可以看出三次样条插值曲线的优越性，即光滑性好，收敛得到保证，局部性好。
3．观察物体的直线运动，得到如下数据
时刻t
0
0.9
1.9
3.0
3.9
5.0
位移s
0
10
30
51
80
111
求运动学方程
。
matlab程序如下：
t=[00.91.93.03.95.0];
s=[010305180111];
n=length（t）;
a2=0;
fori=1:
n
a2=t（i）^2+a2;
end
a3=0;
fori=1:
n
a3=t（i）^3+a3;
end
a4=0;
fori=1:
n
a4=t（i）^4+a4;
end
b1=0;
fori=1:
n
b1=s（i）.*t（i）+b1;
end
b2=0;
fori=1:
n
b2=s（i）.*（t（i）.^2）+b2;
end
A=[a2a3;a3a4]
b=[b1;b2]
c=inv（A）*b
x=0:
0.01:
5.0;
y=c
（1）.*x+c（
升级会员 