大连理工矩阵作业实习题Word下载.docx
《大连理工矩阵作业实习题Word下载.docx》由会员分享，可在线阅读，更多相关《大连理工矩阵作业实习题Word下载.docx（15页珍藏版）》请在冰豆网上搜索。
两种方法均取15位有效数字formatlong，从大到小计算能保证15位有效数字，从小到大只能保证13位。
二、解线性方程组
1．分别利用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组
，其中常向量为
维随机生成的列向量，系数矩阵
具有如下形式
，
其中
为
阶矩阵，
阶单位矩阵，迭代法计算停止的条件为：
，给出n=10,20,30时的不同迭代步数．
构造系数矩阵程序M文件
functionA=xishujuzhen（n）%构造系数矩阵
t=ones（1,n-1）;
t1=ones（1,n-2）;
T=4*diag（t）-diag（t1,1）-diag（t1,-1）;
s=repmat（'
T,'
1,n-1）;
t2=ones（1,（n-1）*（n-2））;
A=eval（sprintf（'
blkdiag（%s）'
s（1:
end-1）））-diag（t2,1-n）-diag（t2,n-1）;
Jacobi法程序M文件
functioni=jacobi（A,x0,b）
i=1;
d=diag（diag（A））;
l=d-tril（A）;
u=d-triu（A）;
d0=inv（d）;
B=d0*（l+u）;
f=d0*b;
x=B*x0+f;
whilenorm（x-x0,inf）>
=1e-10
x0=x;
i=i+1;
i
GS法程序M文件
functioni=GS（a,x0,b）
d=diag（diag（a））;
l=d-tril（a）;
u=d-triu（a）;
B=inv（d-l）*u;
f=inv（d-l）*b;
主程序M文件
formatlong;
n=input（'
n='
）;
b=rand（（n-1）^2,1）;
A=xishujuzhen（n）;
x0=zeros（（n-1）^2,1）;
disp（'
jacobi·
¨
'
i=jacobi（A,x0,b）;
GS·
i=GS（A,x0,b）;
运行结果
n
10
20
30
Jacobi法
i=
431
1738
3909
GS法
220
896
2018
分析
对于此线性方程，Gauss-Seidel迭代法比Jacobi迭代法迭代速度快。
2.用Gauss列主元消去法、QR方法求解如下方程组：
1）Gauss列主元消去法程序
A=[1,2,-1,1,0;
2,5,0,3,4;
1,7,9
2,12;
8,-1,-2,1,-8];
x=zeros（1,4）;
forj=1:
3
k=find（A（j:
end,j）==max（A（j:
end,
j）））+j-1;
a=A（j,:
A（j,:
）=A（k,:
A（k,:
）=a;
fori=（j+1）:
4
A（i,:
）=
）-A（i,j）*A（j,:
）/A（j,j）;
x（1,4）=A（4,5）/A（4,4）;
x（1,3）=（A（3,5）-x（1,4）*A（3,4））/A
（3,3）;
x（1,2）=（A（2,5）-x（1,4）*A（2,4）-x（
1,3）*A（2,3））/A（2,2）;
x（1,1）=（A（1,5）-x（1,4）*A（1,4）-x（
1,3）*A（1,3）-x（1,2）*A（1,2））/A（1,
1）;
x
计算结果
x=（-1,0,1,2）T
2QR程序
B=[1,2,-1,1;
2,5,0,3;
2;
8,-1,-2,1];
c=[0;
4;
12;
-8];
A=B;
Q=eye（4）;
fori=1:
[m,n]=size（A）;
w=A（:
1）-norm（A（:
1））*eye（n,1）;
H=eye（m）-2*w*w'
/（w'
*w）;
a=zeros（m,i-1）;
b=eye（i-1,n+i-1）;
A=H*A;
H=[b;
a,H];
Q=H*Q;
A（:
1）=[];
A（1,:
）=[];
R=Q*B;
d=Q*c;
x（1,4）=d（4,1）/R（4,4）;
x（1,3）=（d（3,1）-x（1,4）*R（3,4））/R
x（1,2）=（d（2,1）-x（1,4）*R（2,4）-x（
1,3）*R（2,3））/R（2,2）;
x（1,1）=（d（1,1）-x（1,4）*R（1,4）-x（
1,3）*R（1,3）-x（1,2）*R（1,2））/R（1,
x=（-1,-1.71766975053489e-015,1,2）T
三、非线性方程的迭代解法
1．求方程
的根，迭代停止的条件为：
；
绘制曲线图
formatlong;
x=0.1:
0.5:
10;
y=exp（x）+2*x.^2+2*sin（x）-log（x）-16;
plot（x,y）;
title（'
曲线图'
由图可知函数单调，在x=2附近有根。
formatlong;
x0=2;
x1=x0-（exp（x0）+2*（x0）^2+2*sin（x0）-log（x0）-16）/（exp（x0）+4*（x0）+2*cos（x0）-1/x）;
whileabs（x1-x0）>
1e-10
x0=x1;
x1
X1=1.96292268319548
选取初值时注意要大于零
2．利用Newton迭代法求多项式
的所有实零点，注意重根的问题。
选取数据点绘制二维曲线图
Matlab程序
x=-3:
y=x.^4-3*x.^3-3*x.^2+11*x-6;
由图像知，在x=-2,x=0.9,x=3处附近存在根，且有一对重根。
Newton法求值程序
x0=input（'
x0='
x1=x0-（x0^4-3*x0^3-3*x0^2+11*x0-6）/（4*x0^3-9*x0^2-6*x0+11）;
0.000001
ifabs（4*x1^3-9*x1^2-6*x1+11）<
0.00001
i=2;
ifabs（12*x1^2-18*x1-6）<
i=3;
x0=-2i=1
x0=0.9i=2
x0=3i=1
利用函数的一阶导、二阶导判断重根的个数。
四、数值积分
用三点Gauss型求积公式计算积分
建立求积公式M文件
functions=guassl（a,b,n）
h=（b-a）/n;
s=0.0;
form=0:
（1*n/2-1）
s=s+h*（guassf（a+h*（（1-1/sqrt（3））+2*m））+guassf（a+h*（（1+1/sqrt（3））+2*m）））;
s;
建立被积函数M文件
functiony=guassf（x）
y=exp（-x）;
运行主程序
s=guassl（0,1,2）
s=
0.631978759531845
五、插值与逼近
1．给定
上的函数
，请做如下的插值逼近：
（1）构造等距节点分别取
的Lagrange插值多项式；
（2）取Chebyshev多项式
的零点：
作插值节点构造
的插值多项式
和上述的插值多项式均要求画出曲线图形（用不同的线型或颜色表示不同的曲线）。
（1）Lagrange程序
t5lag.m
functiony=t5lag（x0,y0,x）
n=length（x0）;
m=length（x）;
fori=1:
m
z=x（i）;
fork=1:
p=1.0;
ifj~=k
p=p*（z-x0（j））/（x0（k）-x0（j））;
s=p*y0（k）+s;
y（i）=s;
主程序
x=linspace（-1,1,1000）;
x5=linspace（-1,1,6）;
x8=linspace（-1,1,9）;
x10=linspace（-1,1,11）;
y=1./（1+x.^2）;
y5=1./（1+x5.^2）;
y8=1./（1+x8.^2）;
y10=1./（1+x10.^2）;
plot（x,y,'
k'
holdon;
yp5=t5lag（x5,y5,x）;
plot（x,yp5,'
b'
yp8=t5lag（x8,y8,x）;
plot（x,yp8,'
g'
yp10=t5lag（x10,y10,x）;
plot（x,yp10,'
r'
legend（'
f（x）'
'
6点插值多项式'
9点插值'
11点插值'
）
n越大插值多项式越接近原函数
（2）主程序
xc=linspace（1,10,10）;
x=cos（（2.*xc-1）./（20）.*pi）;
xp=linspace（-1,0,1000）;
y=1./（1+x.^2）;
yp=t5lag（x,y,xp）;
plot（xp,yp,'
plot（-xp,yp,'
xf=linspace（-1,1,1000）;
yf=1./（1+xf.^2）;
plot（xf,yf,'
Chebyshev'
取Chebyshev多项式
的零点作为插值节点，与原函数更加接近。
3．观察物体的直线运动，得到如下数据
时刻t
0.9
1.9
3.0
3.9
5.0
位移s
51
80
111
求运动学方程
。
clf;
x=[00.91.93.03.95.0];
y=[010305180111];
a=polyfit（x,y,2）;
xd=linspace（0,5,1000）;
yp=a
（1）*xd.^2+a
（2）*xd;
r*'
holdon;
plot（xd,yp）;
original'
polyfit'
二次多项式拟合的比较好。
升级会员 