大连理工大学矩阵与数值分析上机作业13388.docx
《大连理工大学矩阵与数值分析上机作业13388.docx》由会员分享，可在线阅读，更多相关《大连理工大学矩阵与数值分析上机作业13388.docx（26页珍藏版）》请在冰豆网上搜索。
大连理工大学矩阵与数值分析上机作业13388
大连理工大学
矩阵与数值分析上机作业
课程名称：
矩阵与数值分析
研究生姓名：
交作业日时间：
2016年12月20日
第1题
1.1程序：
Clearall;
n=input（'请输入向量的长度n：
'）
fori=1:
n;
v（i）=1/i;
end
Y1=norm（v,1）
Y2=norm（v,2）
Y3=norm（v,inf）
1.2结果
n=10Y1=2.9290
Y2=1.2449
Y3=1
n=100Y1=5.1874
Y2=1.2787
Y3=1
n=1000Y1=7.4855
Y2=1.2822
Y3=1
N=10000Y1=9.7876
Y2=1.2825
Y3=1
1.3分析
一范数逐渐递增，随着n的增加，范数的增加速度减小；二范数随着n的增加，逐渐趋于定值，无群范数都是1.
第2题
2.1程序
clearall;
x
（1）=-10^-15;
dx=10^-18;
L=2*10^3;
fori=1:
L
y1（i）=log（1+x（i））/x（i）;
d=1+x（i）;
ifd==1
y2（i）=1;
else
y2（i）=log（d）/（d-1）;
end
x（i+1）=x（i）+dx;
end
x=x（1:
length（x）-1）;
plot（x,y1,'r'）;
holdon
plot（x,y2）;
2.2结果
2.3分析
红色的曲线代表未考虑题中算法时的情况，如果考虑题中的算法则数值大小始终为1，这主要是由于大数加小数的原因。
第3题
3.1程序
clearall;
A=[1-18144-6722016-40325376-46082304-512];
x=1.95:
0.005:
2.05;
fori=1:
length（x）;
y1（i）=f（A,x（i））;
y2（i）=（x（i）-2）^9;
end
figure（3）;
plot（x,y1）;
holdon;
plot（x,y2,'r'）;
F.m文件
functiony=f（A,x）
y=A
（1）;
fori=2:
length（A）;
y=x*y+A（i）;
end;
3.2结果
第4题
4.1程序
clearall;
n=input（'请输入向量的长度n：
'）
A=2*eye（n）-tril（ones（n,n）,0）;
fori=1:
n
A（i,n）=1;
end
n=length（A）;
U=A;
e=eye（n）;
fori=1:
n-1
[max_data,max_index]=max（abs（U（i:
n,i）））;
e0=eye（n）;
max_index=max_index+i-1;
U=e0*U;
e1=eye（n）;
forj=i+1:
n
e1（j,i）=-U（j,i）/U（i,i）;
end
U=e1*U;
P{i}=e0;%把变换矩阵存到P中
L{i}=e1;
e=e1*e0*e;
end
fork=1:
n-2
Ldot{k}=L{k};
fori=k+1:
n-1
Ldot{k}=P{i}*Ldot{k}*P{i};
end
end
Ldot{n-1}=L{n-1};
LL=eye（n）;
PP=eye（n）;
fori=1:
n-1
PP=P{i}*PP;
LL=Ldot{i}*LL;
end
b=ones（n,2）;
b=e*b;%解方程
x=zeros（n,1）;
x（n）=b（n）/U（n,n）;
fori=n-1:
-1:
1
x（i）=（b（i）-U（i,:
）*x）/U（i,i）;
end
X=U^-1*e^-1*eye（n）;%计算逆矩阵
AN=X';
result2{n-4,1}=AN;
result1{n-4,1}=x;
fprintf（'%d:
\n',n）
fprintf（'%d',AN）;
4.2结果
n=5
1.0625
-0.875
-0.75
-0.5
-0.0625
0.0625
1.125
-0.75
-0.5
-0.0625
0.0625
0.125
1.25
-0.5
-0.0625
0.0625
0.125
0.25
1.5
-0.0625
-0.0625
-0.125
-0.25
-0.5
0.0625
n=10
1.0625
-0.875
-0.75
-0.5
-0.0625
1.0625
-0.875
-0.75
-0.5
-0.0625
0.0625
1.125
-0.75
-0.5
-0.0625
0.0625
1.125
-0.75
-0.5
-0.0625
0.0625
0.125
1.25
-0.5
-0.0625
0.0625
0.125
1.25
-0.5
-0.0625
0.0625
0.125
0.25
1.5
-0.0625
0.0625
0.125
0.25
1.5
-0.0625
-0.0625
-0.125
-0.25
-0.5
0.0625
-0.0625
-0.125
-0.25
-0.5
0.0625
1.0625
-0.875
-0.75
-0.5
-0.0625
1.0625
-0.875
-0.75
-0.5
-0.0625
0.0625
1.125
-0.75
-0.5
-0.0625
0.0625
1.125
-0.75
-0.5
-0.0625
0.0625
0.125
1.25
-0.5
-0.0625
0.0625
0.125
1.25
-0.5
-0.0625
0.0625
0.125
0.25
1.5
-0.0625
0.0625
0.125
0.25
1.5
-0.0625
-0.0625
-0.125
-0.25
-0.5
0.0625
-0.0625
-0.125
-0.25
-0.5
0.0625
同样的方法可以算出n=20,n=30时的结果，这里就不罗列了。
第5题
5.1程序
clearall;
n=input（'请输入向量的长度n：
10至20'）
fori=1:
n
forj=1:
n
a（i,j）=1/（i+j-1）;
end
end
forj=1:
n
sum=0;
fork=1:
j-1
sum=sum+l（j,k）^2;
end
l（j,j）=sqrt（a（j,j）-sum）;
fori=j+1:
n
sum=0;
fork=1:
j-1
sum=sum+l（i,k）*l（j,k）;
end
l（i,j）=（a（i,j）-sum）/l（j,j）;
end
end
b=ones（n,1）;
y=zeros（n,1）;
y（n）=b（n）/l（n,n）;
fori=n-1:
-1:
1
y（i）=（b（i）-l（i,:
）*y）/l（i,i）;
end
l=l';
x=zeros（n,1）;
x（n）=y（n）/l（n,n）;
fori=n-1:
-1:
1
x（i）=（y（i）-l（i,:
）*x）/l（i,i）;
end
fprintf（'%d\t',x）;
fprintf（'\n'）;
5.2结果
n=10
n=11
n=12
n=13
n=14
n=15
n=16
n=17
n=18
n=19
n=20
-746517.8688
3111493.423
-11884558.85
478355909.6
497329749.5
519377549.9
445748685
378885969.7
341571897.8
1020943960
994971030.6
69820595.08
-354658479.9
1634110905
-77484610111
-80547115920
-82914903907
-72327488690
-62111481037
-60010758946
-1.76915E+11
-1.68795E+11
-1587444197
9875543407
-54993461592
3.06265E+12
3.18327E+12
3.23484E+12
2.86169E+12
2.48062E+12
2.55037E+12
7.40163E+12
6.93865E+12
152********
-1.17236E+11
7.93546E+11
-5.20373E+13
-5.40791E+13
-5.42806E+13
-4.8713E+13
-4.26882E+13
-4.65382E+13
-1.3123E+14
-1.21156E+14
-76184620048
7.35352E+11
-6.11155E+12
4.7524E+14
4.93812E+14
4.89563E+14
4.46792E+14
3.97243E+14
4.57473E+14
1.22E+15
1.11E+15
2.18036E+11
-2.70378E+12
2.80302E+13
-2.62E+15
-2.72E+15
-2.66E+15
-2.48E+15
-2.25E+15
-2.72E+15
-6.64E+15
-5.98E+15
-3.70513E+11
6.12295E+12
-8.11E+13
9.24E+15
9.59E+15
9.27E+15
8.90E+15
8.30E+15
1.04E+16
2.17E+16
1.94E+16
3.69292E+11
-8.64269E+12
1.51789E+14
-2.17E+16
-2.25E+16
-2.14E+16
-2.14E+16
-2.07E+16
-2.60E+16
-4.18E+16
-3.74E+16
-1.99261E+11
7.40507E+12
-1.83339E+14
3.40E+16
3.53E+16
3.34E+16
3.52E+16
3.56E+16
4.28E+16
4.06E+16
3.77E+16
44905979430
-3.52275E+12
1.3792E+14
-3.55E+16
-3.68E+16
-3.47E+16
-3.89E+16
-4.12E+16
-4.38E+16
-4.89E+15
-1.01E+16
7.13565E+11
-5.87483E+13
2.35E+16
2.44E+16
2.37E+16
2.73E+16
2.94E+16
2.40E+16
-1.64E+16
-5.28E+15
1.08203E+13
-8.98E+15
-9.29E+15
-1.03E+16
-9.10E+15
-7.78E+15
-4.50E+15
-1.23E+16
-1.58E+16
1.50E+15
1.55E+15
3.06E+15
-2.94E+15
-6.68E+15
1.72E+15
1.63E+16
1.04E+16
2.17168E+12
-7.94676E+14
5.08E+15
7.29E+15
-1.45E+15
3.84E+16
3.85E+16
1.58892E+14
-2.50E+15
-2.19E+15
-8.57E+15
-2.64E+16
-2.01E+16
4.80384E+14
-3.1296E+14
1.48E+16
-8.62E+16
-7.98E+16
2.304E+14
-9.01E+15
1.40E+17
1.19E+17
1.99E+15
-8.07E+16
-6.35E+16
1.70E+16
1.09E+16
7.5453E+14
第6题
6.1程序
clearall;
A=[1234;
-13sqrt
（2）sqrt（3）;
-22exp
（1）pi;
-sqrt（10）2-37;
0275/2;];
U=f61（A（:
2））;
HA=f62（U,A）;
f.m文件
functionU=f61（x）
e1=eye（length（x）,1）;
U=x-sign（x
（1））*sqrt（dot（x,x））*e1;
U=U./sqrt（dot（U,U））;
functionHA=f62（U,A）
HA=A-2*U*U'*A;
6.2结果
-2.264911064
5
4.735840869
7.695867546
2.264911064
4.44E-16
-0.321627306
-1.963816738
0.176607376
4.44E-16
1.561054583
0.677680956
-0.985670284
4.44E-16
-4.157227246
4.536088303
2.176607376
4.44E-16
5.842772754
0.036088303
第7题
7.1程序
clearall;
max=1000;
x（1,:
）=[123];
fori=1:
max
x（i+1,1）=（-2+x（i,2）+3*x（i,3））/5;
x（i+1,2）=（1+x（i,1）-4*x（i,3））/2;
x（i+1,3）=（10+3*x（i,1）-4*x（i,2））/15;
err（i）=sqrt（dot（x（i+1,:
）-x（i,:
）,x（i+1,:
）-x（i,:
）））;
if（err（i）<10^-6）
break;
end
end
figure（7）;
plot（err）;
clearerr;
x（1,:
）=[123];
fori=1:
max
x（i+1,1）=（-2+x（i,2）+3*x（i,3））/5;
x（i+1,2）=（1+x（i+1,1）-4*x（i,3））/2;
x（i+1,3）=（10+3*x（i+1,1）-4*x（i+1,2））/15;
err（i）=sqrt（dot（x（i+1,:
）-x（i,:
）,x（i+1,:
）-x（i,:
）））;
if（err（i）<10^-6）
break;
end
end
holdon
plot（err,'r'）;
7.2结果
误差越来越小。
第8题
8.1程序
clearall;
max=100;
x
（1）=1;
fori=1:
max
x（i+1）=x（i）-（x（i）^3+2*x（i）^2+10*x（i）-100）/（3*x（i）^2+4*x（i）+10）;
if（abs（x（i+1）-x（i））<10^-6）
break;
end
end
figure（8）
plot（x）;
clearx;
x
（1）=0;
x
（2）=1;
fori=2:
max
x（i+1）=x（i）-（x（i）^3+2*x（i）^2+10*x（i）-100）/（x（i）^3+2*x（i）^2+10*x（i）-（x（i-1）^3+2*x（i-1）^2+10*x（i-1）））*（x（i）-x（i-1））;
if（abs（x（i+1）-x（i））<10^-6）
break;
end
end
holdon
plot（x,'r'）;
8.2结果
8.3分析
由计算结果可知，弦截法的收敛速度比牛顿法的收敛速度快。
第9题
9.1程序
clearall;
f（0,4*pi）;
f.m文件
function[]=f（l,r）
if（r-l<10^-6）
fprintf（'%g,',（r-l）/2+l）;
return
end
iff9（l）*f9（r）<0
iff9（（l+r）/2+l）*f9（l）<0
f（l,（l+r）/2+l）;
end
iff9（（l+r）/2+l）*f9（r）<0
f（（l+r）/2+l,r）;
end
end
iff9（l）*f9（r）>0
f（l,（l+r）/2+l）;
f（（l+r）/2+l,r）;
end
9.2结果
X=4.71239、8.24668、17.6715、14.1372
第10题
10.1程序
clearall;
n=3;%节点个数
Xj=0:
1/n:
1;
y=sin（pi*Xj）;
fori=1:
n+1
f（i,1）=y（i）;
end
forj=2:
n+1
fori=1:
n-j+2;
dx=（（j-1）/n）;
f（i,j）=（f（i+1,j-1）-f（i,j-1））/dx;
end
end
fori=1:
n+1
a（i）=f（1,i）;
end
x=0:
0.001:
1;
fori=1:
1/0.001+1;
y1（i）=sin（pi*x（i））;
y2（i）=f10（a,Xj,x（i））;
end
figure（10）;
plot（x,y1,'r'）;
holdon;
plot（x,y2,''）;
10.2结果
10.3分析
有图像可知插值函数的值已经很接近原函数的值了。
第11题
11.1程序
clearall;
n=input（'请输入n：
'）%n代表节点
Xj=-5:
1/n:
5;
Yj=1./（1.+Xj.^2）;
x=-5:
0.01:
5;
fori=1:
10/0.01+1;
y1（i）=1/（1+x（i）^2）;
y2（i）=f（Yj,Xj,x（i））;
end
figure（11）;
plot（x,y1,'r'）;
holdon;
plot（x,y2,,''）;
f.m文件
functiony=f（Yj,Xj,x）
y=0;
fori=1:
length（Yj）
l=1;
forj=1:
length（Xj）
ifi==j
continue;
end
l=l*（x-Xj（j））/（Xj（i）-Xj（j））;
end
y=y+Yj（i）*l;
end
11.2结果
从左往右n依次取1、2、3、4
11.3分析
随着n的不断增加，插值越来越接近真实值。
第12题
12.1程序
clearall;
n=input（'请输入n：
'）%n=501002005001000
x=0:
2*pi/n:
2*pi;
fori=1:
n+1
X=x（i）;
Y（i）=exp（3*X）*cos（pi*X）;
end
Y
（1）=1;
fori=1:
n
X=（x（i）+x（i+1））/2;
Y2（i）=exp（3*X）*cos（pi*X）;
end
T=（x（n+1）-x
（1））/（2*n）*（Y
（1）+2*sum（Y（2:
n））+Y（n+1））;
S=（x（n+1）-x
（1））/（6*n）*（Y
（1）+4*sum（Y2）+2*sum（Y（2:
n））+Y（n+1））;
13.2结果
n
50
100
200
500
1000
T
3.512534119493697e+07
3.520489119998542e+07
3.522553497836321e+07
3.523136936592422e+07
3.523220478232005e+07
S
3.523140786833490e+07
3.523241623782247e+07
3.523247916807637e+07
3.523248325445201e+07
3.523248335509114e+07
第13题
13.1程序
clearall;
G2=f（1/sqrt（3））+f（-1/sqrt（3））;
G3=0.555555556*f（-0.7745966692）+0.555555556*f（+0.7745966692）+0.88888888889*f（0）;
G5=0.2369268851*f（-0.9061798459）+0.2369268851*f（+0.9061798459）+...0.4786286705*f（-0.5384693101）+0.4786286705*f（+0.5384693101）+...0.5688888889*f（0）;
f.m文件1
functiony=f（x）
y=x^2/sqrt（1-x^2）;
f.m文件2
functiony=f（x）
x=pi/2*（x+1）/2;
y=sin（x）/x;
13.2结果
（1）y=x^2/sqrt（1-x^2）
（2）y=sin（x）/x
G2
0.816496*********
G2
1.74487172393660
G3
1.05409255403515
G3
1.74531021550858
G5
1.05409255403515
G5
1.74530859901330
第14题
14.1程序
clearall;
U0=2;%初值
step=0.2;%step表示步长
U1
（1）=U0;
Len=1;
fori=1:
floor（Len/step）
U1（i+1）=U1（i）+step*f14（U1（i）,（i-1）*step）;
end
%改进的Euler法
U2
（1）=U0;
fori=1:
floor（Len/step）
U2（i+1）=U2（i）+step/2*（f14（U2（i）,（i-1）*step）+f14（U1（i+1）,（i+1-1）*step））;
end
%Runge_KUtta
U3
（1）=U0;
fori=1:
floor（Len/step）
t=（i-1）*ste
升级会员 