大连理工矩阵作业实习题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'
二次多项式拟合的比较好。