2016矩阵与数值分析上机作业满分.docx
《2016矩阵与数值分析上机作业满分.docx》由会员分享，可在线阅读，更多相关《2016矩阵与数值分析上机作业满分.docx（25页珍藏版）》请在冰点文库上搜索。
第一题
考虑计算给定向量的范数：
输入向量x=x1,x2,…,xnT,输出x1,x2,x∞。
请编制一个通用程序，并用你编制的程序计算如下向量的范数：
x=1,12,13,…,1nT,y=1,2,…,nT.
对n=10，100，1000甚至更大的n计算其范数，你会发现什么结果?
你能否修改你的程序使得计算结果相对精确呢？
代码：
functionProblem1（）
clc
N=input（'inputN='）;
X=zeros（1,N）;
Y=1:
N;
fori=1:
N
X（i）=1/i;
end
x1=0;
x2=0;
y1=0;
y2=0;
fori=1:
N
x1=x1+abs（X（i））;
x2=x2+X（i）*X（i）;
y1=y1+abs（Y（i））;
y2=y2+Y（i）*Y（i）;
end
x1;
x2=sqrt（x2）;
x3=max（abs（X））;
y1;
y2=sqrt（y2）;
y3=max（abs（Y））;
fprintf（'x的1范数是%.2f\n',x1）;
fprintf（'x的2范数是%.2f\n',x2）;
fprintf（'x的无穷范数是%.2f\n',x3）;
fprintf（'y的1范数是%.2f\n',y1）;
fprintf（'y的2范数是%.2f\n',y2）;
fprintf（'y的无穷范数是%.2f\n',y3）;
结果：
inputN=10
x的1范数是2.93
x的2范数是1.24
x的无穷范数是1.00
y的1范数是55.00
y的2范数是19.62
y的无穷范数是10.00
inputN=100
x的1范数是5.19
x的2范数是1.28
x的无穷范数是1.00
y的1范数是5050.00
y的2范数是581.68
y的无穷范数是100.00
inputN=1000
x的1范数是7.49
x的2范数是1.28
x的无穷范数是1.00
y的1范数是500500.00
y的2范数是18271.11
y的无穷范数是1000.00
第二题
考虑y=fx=ln（1+x）x,其中定义f（0）=1，此时f（x）是连续函数。
用此公式计算当x∈-10-15,10-15时的函数值，画出图像。
另一方面，考虑下面算法：
d=1+x
ifd=1then
y=1
else
y=lnd/（d-1）
endif
用此算法计算x∈-10-15,10-15时的函数值，画出图像。
比较一下发生了什么？
代码：
functionproblem2（）
clc
N=1000;
n=2*10^（-15）/N;
%公式计算
f=zeros（1,N+1）;
t=1;
fori=-10^（-15）:
n:
10^（-15）
if（i~=0）
f（t）=log（1+i）/i;
else
f（t）=1;
end
t=t+1;
end
%算法计算
a=zeros（1,N+1）;
t=1;
fori=-10^（-15）:
n:
10^（-15）
d=1+i;
if（d~=1）
a（t）=log（d）/（d-1）;
else
a（t）=1;
end
t=t+1;
end
%画图，左边是公式算的，右边是算法算的
subplot（1,2,1）;
plot（-10^（-15）:
n:
10^（-15）,f）;
holdon
subplot（1,2,2）;
plot（-10^（-15）:
n:
10^（-15）,a）;
holdon
结果：
第三题
首先编制一个利用秦九韶算法计算一个多项式在给定点的函数值的通用程序，你的程序包括输入多项式的系数以及给定点，输出函数值。
利用你编制的程序计算
px=x-29=x9-18x8+144x7-672x6+2016x5-4032x4+5376x3-4608x2-512在x=2邻域附近的值。
画出px在x∈1.95,20.5上的图像。
代码：
functionproblem3（）
%秦九韶算法
clc
n=9;
c=[1,-18,144,-672,2016,-4032,5376,-4608,2304,-512];
N=100;
m=（2.05-1.95）/N;
y=zeros（1,N+1）;
ydex=0;
forx=1.95:
m:
2.05
ydex=ydex+1;
dex=1;
temp=c（dex）;
whiledex<=n
dex=dex+1;
temp=temp*x+c（dex）;
end
y（ydex）=temp;
end
plot（1.95:
m:
2.05,y）;
结果：
第四题
编制计算给定矩阵A的LU分解和PLU分解的通用程序，然后用你编制的程序完成下面两个计算任务：
（1）考虑
A=10⋯-1⋱⋱⋮⋱⋱0⋮01⋮1-1⋯-111-1⋯-1-11∈Rn×n,
自己取定x∈Rn，并计算b=Ax。
然后用你编制的不选主元和列主元的Gauss消去法求解该方程组，记你计算出的解为x。
对n从5到30估计计算解的精度。
（2）对n从5到30计算其逆矩阵。
代码：
%LU分解
functionproblem4（）
clc
%生成矩阵A
forn=5:
30
A=-ones（n）;
A（:
n）=ones（n,1）;
foria=1:
n
A（ia,ia）=1;
forja=（ia+1）:
（n-1）
A（ia,ja）=0;
end
end
x=rand（n,1）;
b=A*x;
%LU分解
[n,n]=size（A）;
L=eye（n）;
M=eye（n）;
fori=1:
n-1
ifA（i,i）==0
fprintf（'Error:
A（%d,%d）=0!
\n',i,i）;return;
end
forj=i+1:
n
L（j,i）=A（j,i）/A（i,i）;
end
fork=i+1:
n
forl=i:
n
A（k,l）=A（k,l）-L（k,i）*A（i,l）;
end
end
end
U=A;
M=inv（U）*inv（L）*M
y1=inv（L）*b;
x1=inv（U）*y1;
error=norm（x-x1）
end
return;
%PLU分解
functionproblem4（）
clc
%生成矩阵A
n=5;
A=-ones（n）;
A（:
n）=ones（n,1）;
foria=1:
n
A（ia,ia）=1;
forja=（ia+1）:
（n-1）
A（ia,ja）=0;
end
end
x=rand（n,1）;
b=A*x;
%PLU分解
[n,n]=size（A）;
Q=eye（n）;
R=eye（n）;
S=eye（n）;
L=zeros（n,n,n-1）;
foril=1:
（n-1）
L（:
:
il）=eye（n）;
end
P=zeros（n,n,n-1）;
foril=1:
（n-1）
P（:
:
il）=eye（n）;
end
fori=1:
n-1
m=i;
a=abs（A（i,i））;
fork=i+1:
n
ifabs（A（k,i））>a
a=abs（A（k,i））;
m=k;
end
end
A（[i,m],:
）=A（[m,i],:
）;
P（[i,m],:
i）=P（[m,i],:
i）;
forj=i+1:
n
L（j,i,i）=-A（j,i）/A（i,i）;
forl=i:
n
A（j,l）=A（j,l）+L（j,i,i）*A（i,l）;
end
end
end
U=A;
U
forx=1:
n-1
Q=L（:
:
x）;
fory=x+1:
n-1
Q=P（:
:
y）*Q*P（:
:
y）;
end
S=S*inv（Q）;
end
%L矩阵
S
forz=1:
（n-1）
R=P（:
:
z）*R;
end
%P矩阵
R
b1=R*b;
b2=inv（S）*b1;
x2=inv（U）*b2;
M=inv（U）*inv（S）*R
error=norm（x-x2）
结果：
由于篇幅有限仅列出n=5时的精度及逆矩阵
LU分解
M=
0.5000-0.2500-0.1250-0.0625-0.0625
00.5000-0.2500-0.1250-0.1250
000.5000-0.2500-0.2500
0000.5000-0.5000
0.50000.25000.12500.06250.0625
error=
8.8818e-16
PLU分解
M=
0.5000-0.2500-0.1250-0.0625-0.0625
00.5000-0.2500-0.1250-0.1250
000.5000-0.2500-0.2500
0000.5000-0.5000
0.50000.25000.12500.06250.0625
error=
7.8516
第五题
编制计算对称正定阵的Cholesky分解的通用程序，并用你编制的程序计算Ax=b，其中A=aijϵRn×n,aij=1i+j-1。
b可以由你自己取定，对n从10到20验证程序的可靠性。
代码：
functionproblem5（）
%cholesky分解
clc
%生成矩阵A
N=10;
A=zeros（N）;
foria=1:
N
forja=1:
N
A（ia,ja）=1/（ia+ja-1）;
end
end
L=chol（A）;
L
结果：
由于篇幅有限仅列出n=10时的L矩阵
L=
1.00000.50000.33330.25000.20000.16670.14290.12500.11110.1000
00.28870.28870.25980.23090.20620.18560.16840.15400.1417
000.07450.11180.12780.13310.13310.13040.12650.1220
0000.01890.03780.05250.06300.07020.07480.0777
00000.00480.01190.01950.02650.03260.0378
000000.00120.00360.00680.01030.0139
0000000.00030.00110.00220.0038
00000000.00010.00030.0007
000000000.00000.0001
0000000000.0000
第六题
（1）编制程序House（x），其作用是对输入的向量x，输出单位向量u使得I-2uuTx=x2e1。
（2）编制Householder变换阵H=I-2uuTϵRn×n乘以AϵRn×m的程序HA，注意，你的程序并不显式的计算出H。
（3）考虑矩阵
A=123-132-22e43π-102-302775/2,
用你编制的程序计算H使得HA的第一列为αe1的形式，并将HA的结果显示。
代码：
functionproblem6（）
clc
x=[1;-1;-2;-10^（1/2）;0];
ifx
（1）>0
sg=1;
else
sg=-1;
end
Nx=length（x）;
e1=zeros（Nx,1）;
e1
（1）=1;
sum=0;
fori=1:
Nx
sum=sum+x（i）*x（i）;
end
x2=sqrt（sum）;
clearsum;
w=x-sg*x2*e1;
u=w/sqrt（w'*w）;
u
%生成矩阵A
A=[1,2,3,4;-1,3,2^（1/2）,3^（1/2）;-2,-2,exp
（1）,pi;-10^（1/2）,2,-3,7;0,2,7,5/2];
H=eye（Nx）-2*u*u';
H*A
结果：
u=
-0.6124
-0.2041
-0.4082
-0.6455
0
ans=
4.0000-0.83111.4090-6.5378
0.00002.05630.8839-1.7805
0.0000-3.88741.6576-3.8836
0.0000-0.9843-4.6770-4.1078
02.00007.00002.5000
第七题
用Jacobi和Gauss-Seidel迭代求解下面的方程组，输出迭代每一步的误差xk-x∅:
5x1-x2-3x3=-2-x1+2x2+4x3=1-3x1+4x2+15x3=10
代码：
functionproblem7（）
clc
inter=10;%迭代次数
rx=[69*7/61-8,69*3/（61*2）-7/2,69/61]';
%根据题意生成A,b
A=[5-1-3;-124;-3415];
b=[-2110]';
N=length（b）;
x=ones（N,1）;
%Jacobi迭代法
L=zeros（N,N）;
D=L;
U=L;
foria=1:
N
D（ia,ia）=A（ia,ia）;
ifiaforiu=（ia+1）:
N
U（ia,iu）=A（ia,iu）;
end
end
ifia>1
foril=1:
（ia-1）
L（ia,il）=-A（ia,il）;
end
end
end
forintt=1:
inter
x=inv（D）*（L+U）*x+inv（D）*b;
wuca1=0;
%求误差
forix=1:
N
wuca1=wuca1+（x（ix）-rx（ix））^2;
end
wuca1
end
%G-S迭代法
forintt=1:
inter
x=inv（D-L）*U*x+inv（D-L）*b;
wuca2=0;
%求误差
forix=1:
N
wuca2=wuca2+（x（ix）-rx（ix））^2;
end
wuca2
end
结果:
wuca1=
24.6036
wuca1=
12.3259
wuca1=
1.9276
wuca1=
5.4520
wuca1=
15.9787
wuca1=
10.0748
wuca1=
3.6731
wuca1=
6.2770
wuca1=
11.9670
wuca1=
8.9591
wuca2=
5.5096
wuca2=
9.1425
wuca2=
6.9372
wuca2=
8.0757
wuca2=
7.4675
wuca2=
7.7807
wuca2=
7.6170
wuca2=
7.7018
wuca2=
7.6577
wuca2=
7.6806
第八题
取不同的初值用Newton迭代以及弦截法求方程x3+2x2+10x-100=0的实根，列表或者画图说明收敛速度。
代码：
functionproblem8（）
%inNew记录Newton迭代法的收敛速度
%inS记录弦截法的收敛次数
clc
chuzhi=4;%初值
f=[];
%Newton迭代
x=chuzhi;
forinNew=1:
inter
xx=x-（x^3+2*x^2+10*x-100）/（3*x^2+4*x+10）;
f（inNew）=x;
ifabs（xx-x）<10^（-5）
break;
end
x=xx;
end
subplot（1,2,1）;
plot（1:
inNew,f）;
holdon
clearf;
inNew
xx
%弦截法
xk=chuzhi;
xkr=chuzhi-1;
forinS=1:
inter
xka=xk-（xk^3+2*xk^2+10*xk-100）*（xk-xkr）/（（xk^3+2*xk^2+10*xk-100）-（xkr^3+2*xkr^2+10*xkr-100））;
f（inS）=xk;
ifabs（xka-xk）<10^（-5）%精度10^（-5）
break;
end
xkr=xk;
xk=xka;
end
subplot（1,2,2）;
plot（1:
inS,f）;
holdon
inS
xka
结果：
inNew=
7
xx=
3.4606
inS=
10
xka=
3.4606
结果：
第九题
用二分法求方程excosx+2=0在区间0,4π上的所有根。
代码：
functionproblem9（）
%二分法求解方程
clc
%查看根的个数
N=50;
yk=zeros（1,N+1）;
h=4*pi/N;
e=[];
c=[];
i=1;
forti=0:
h:
4*pi
e（i）=-2/（exp（ti））;
c（i）=cos（ti）;
i=i+1;
end
plot（0:
h:
4*pi,e）;
holdon
plot（0:
h:
4*pi,c）;
%根的取值范围
min=[14710];
max=[25811];
accr=zeros（1,4）;
acc=0.0001;%精度
forr=1:
4
while（max（r）-min（r）>acc）
ymax=exp（max（r））*cos（max（r））+2;
temp=（max（r）+min（r））/2;
ytemp=exp（temp）*cos（temp）+2;
ifytemp*ymax==0
break;
else
ifytemp*ymax>0
max（r）=temp;
else
min（r）=temp;
end
end
end
accr（r）=temp;
end
accr
结果：
accr=
1.88074.69407.854810.9955
第十题
考虑函数f（x）=sin（πx），x∈[0,1]。
用等距节点作f（x）的Newton插值，画出插值多项式以及f（x）的图像，观察收敛性。
代码：
functionf=Newton（x,y）
symst;
n=length（x）;
c（1:
n）=0;
f=y
（1）;
y1=0;
p=1;
fori=1:
n-1
forj=i+1:
n
y1（j）=（y（j）-y（i））/（x（j）-x（i））;
end
c（i）=y1（i+1）;
p=p*（t-x（i））;
f=f+c（i）*p;
y=y1;
end
f=simplif
升级会员 