2016级矩阵与数值分析上机作业.docx
《2016级矩阵与数值分析上机作业.docx》由会员分享，可在线阅读，更多相关《2016级矩阵与数值分析上机作业.docx（72页珍藏版）》请在冰点文库上搜索。
2016级矩阵与数值分析上机作业
学生班级：
学生姓名：
任课教师：
所在学院：
电子信息与电气学部学生学号：
使用软件：
MATLAB
2016年12月10号
1.考虑计算给定向量的范数：
输入向量x=（x,x,...,x）T，输出x，x，
1 2 n 1 2
x¥请编制一个通用程序，并用你编制的程序计算如下向量的范数：
ç
x=æ1,
è
1,1
23
,...,
1öT
n
÷ ，
ø
y=（1,2,3,...,n）T
对n=10，100，1000甚至更大的n计算其范数，你会发现什么结果？
你能否修改你的程序使得计算结果相对精确呢？
（1）计算范数的程序
function[Nor1,Nor2,Nor3]=Nor（a）b=length（a）;
formatlongNor11=0;
formatlongNor21=0;i=1;
whilei<=b
Nor11=Nor11+abs（a（i））;Nor21=Nor21+a（i）^2;t=a
（1）;
ifabs（a（i））>t
t=abs（a（i））;
endi=i+1;end
Nor1=Nor11; %计算1范数
Nor2=sqrt（Nor21）; %计算2范数Nor3=t; %计算无穷范数end
（2）x与y向量的程序
functiona=afun（n）a=zeros（n）;
a
（1）=1;
fori=1:
1:
n;a（i）=1/i;
endend
functionb=bfun（n）b=zeros（n）;
b
（1）=1;
fori=1:
1:
n;b（i）=i;
endend
运行：
>>
n=10;
>>
a=afun（n）;
>>
b=bfun（n）;
>>
[f1,f2,f3]=Nor（a）;
>>
[h1,h2,h3]=Nor（b）;
f1
=2.928968253968254 f2=1.244896674895769
f3=
1
h1
=55 h2=19.621416870348583
h3=
10
>>n=100;
a=afun（n）;b=bfun（n）;[f1,f2,f3]=Nor（a）;
>>[h1,h2,h3]=Nor（b）;
f1=5.187377517639621 f2=1.278664889713052 f3=1h1=5050 h2=5.816786054171153e+02 h3=100
>>n=1000;
a=afun（n）;b=bfun（n）;
[f1,f2,f3]=Nor（a）;
>>[h1,h2,h3]=Nor（b）;
f1=7.485470860550343 f2=1.282160117411847 f3=1h1=500500 h2=1.827111107732642e+04 h3=1000
上述结果由于MATLAB的高精度运算较为准确，但是由于当n非常大时可能会出现了大
数吃小数的现象使误差增大，而y向量的范数结果较为准确是由于y向量是按从小到大排列。
改进范数计算程序使输入序列按从小到大排列再计算其范数：
function[Nor1,Nor2,Nor3]=Nor（c）a=sort（c）;b=length（a）;
formatlongNor11=0;
formatlongNor21=0;i=1;
whilei<=b
Nor11=Nor11+abs（a（i））;Nor21=Nor21+a（i）^2;t=a
（1）;ifabs（a（i））>t
t=abs（a（i））;
endi=i+1;end
Nor1=Nor11; %计算1范数
Nor2=sqrt（Nor21）; %计算2范数
Nor3=t; %计算无穷范数
end
2.考虑y=
f（x）=ln（x+1），其中定义f（0）=1，此时f（x）是连续函数。
用此公式计算当
x
xÎ[-10-15,10-15
时的函数值，画出图像。
另一方面，考虑下面算法：
d=x+1ifd=1theny=1
else
y=lnd/（d-1）
endif
用此算法计算x∈[−10−15,10−15]时的函数值，画出图像。
比较一下发生了什么？
函数：
functionF=Piecewise1_x（x）
F=log（x+1）/x*（x<0）+1*（x>=0&x<=0）+log（x+1）/x*（x>0）;end
functionF=Piecewise2_x（x）
F=log（x）/（x-1）*（x<1）+1*（x>=1&x<=1）+log（x）/（x-1）*（x>1）;end
运行：
clcclear
x=linspace（-10^（-15）,10^（-15））;d=x+1;
%原来图像（红色曲线）
F=Piecewise1_x（x）;%计算相应函数值plot（x,F,'r'）;%绘制曲线
holdon;
%改变算法后的图像（蓝色曲线）
F1=Piecewise2_x（d）;%计算相应函数值plot（x,F1,'b'）;%绘制曲线
holdon;
plot（0*ones（1,2）,ylim,'g:
'）;%画区间间隔线
plot（xlim,1*ones（1,2）,'g:
'）;%画区间间隔线xlabel（'变量X'）
ylabel（'变量Y1&Y2'）
由上图可知改变后的算法画出的曲线（蓝色）比直接画出的曲线（红色）更贴近实际值误差更小，且由于ln（x+1）<=x所以不会出现小数做除数导致误差增大。
这可能是因为在原算法中x+1出现了大数吃小数的现象导致的结果。
3.首先编制一个利用秦九韶算法计算一个多项式在给定点的函数值的通用程序，
你的程序包括输入多项式的系数以及给定点，输出函数值。
利用你编制的程序计算
p（x）=（x-2）9=x9-18x8+144x7-672x6+2016x5-4032x4+5376x3-4608x2+2304x-512
在x=2邻域附近的值。
画出p（x）在x∈[1.95,20.5]上的图像。
秦九韶算法程序：
function p=p（a,x）n=length（a）;p=a（n）;
fori=n:
-1:
2pi=p.*x+a（i-1）;p=pi;
endend
运行：
clcclear
x=linspace（1.95,20.5）;
a=[-512,2304,-4608,5376,-4032,2016,-672,144,-18,1];
F=p（a,x）;%计算相应函数值plot（x,F,'r'）;%绘制曲线holdon;
4.编制计算给定矩阵A的LU分解和PLU分解的通用程序，然后用你编制的程序完成下面两个计算任务：
（1）考虑：
é1 0
ê
ê-1O
A=êM O
ê
ê-1L
êë-1L
L 0
O M
O 0
-1 1
-1 -1
1ù
M
ú
ú
1ú
ú
1ú
1úû
自己取定xÎRn，并计算b=Ax，然后利用你编制的不选主元和列主元的
L
Gausss消去法求解该方程组，记你算的解为x，对n从5到30估计计算解得精
度。
（2）对n从5到30计算其逆矩阵
（1）
1.LU分解
%LU分解通用程序
function[L,U]=zlu（A）
%ZLU-LUdecompositionformatrixA
%workasgausselimination[m,n]=size（A）;
ifm~=n
error（'Error,currenttimeonlysupportsquarematrix'）;
end
L=zeros（n）;
U=zeros（n）;
fork=1:
n
gauss_vector=A（:
k）;
gauss_vector（k+1:
end）=gauss_vector（k+1:
end）./gauss_vector（k）;
gauss_vector（1:
k）=zeros（k,1）;
L（:
k）=gauss_vector;
L（k,k）=1;
forl=k+1:
n
A（l,:
）=A（l,:
）-gauss_vector（l）*A（k,:
）;
end
end
U=A;
%LU误差分析程序
function d=zlu1（n）x1=zeros（n,1）;x=zeros（n,1）;
fori=1:
nx1（i）=i;
endA=zeros（n,n）;fori=1:
1:
n
forj=1:
1:
nifi>jA（i,j）=-1;
elseif （ielseA（i,j）=1;
endend
endb=A*x1;
[L,U]=zlu（A）;
y=zeros（n,1）;
y
（1）=b
（1）;
fori=2:
n
y（i）=b（i）-sum（L（i,1:
i-1）'.*y（1:
i-1））;end
%yx（n）=y（n）/U（n,n）;fori=n-1:
-1:
1
x（i）=（y（i）-sum（U（i,i+1:
n）'.*x（i+1:
n）））/U（i,i）;end
%x
d=norm（x1-x）;%计算算法误差end
运行：
cleard=zeros（26,1）;forn=5:
30
d（n-4）=zlu1（n）;
endd1=d';d1
d1=
Columns1through9
0 0 0 0 0 0 0 0 0
Columns10through18
0 0 0 0 0 0 0 0 0
Columns19through26
0 0 0 0 0 0 0 0
由程序计算结果可知，LU分解计算程序计算误差非常小，基本无误差计算
2.PLU分解
%PLU分解,PA=LU
function[P,L,U]=PLU（A）formatratB=A;[m,n]=size（B）;
Q=[1:
m]';P=zeros（m）;x=zeros（1,m）;
fori=1:
m
forj=i+1:
m
%--------交换行------
ifB（i,i）==0fork=i+1:
m
ifB（k,i）~=0x=B（k,:
）;B（k,:
）=B（i,:
）;B（i,:
）=x;t=Q（k）;Q（k）=Q（i）;Q（i）=t;
break
end
end
end
%--------交换行------
fork=m:
-1:
i+1
B（j,k）=B（j,k）-B（j,i）/B（i,i）*B（i,k）;
endB（j,i）=B（j,i）/B（i,i）;
end
end
fori=1:
m
s=Q（i）;
P（i,s）=1;
end
L=eye（m）+tril（B,-1）;U=triu（B）;
%fprintf（'PLU分解矩阵：
\n以下为结果：
\n'）
%Q,B,P,L,U
%fprintf（'原矩阵：
A='）;pretty（sym（A））;
%fprintf（'变换系数：
Q='）;pretty（sym（Q））;
%fprintf（'变换矩阵：
P='）;pretty（sym（P））;
%fprintf（'下三角阵：
L='）;pretty（sym（L））;
%fprintf（'上三角阵：
U='）;pretty（sym（U））;
%LU误差分析程序
functiond=PLU1（n）x1=zeros（n,1）;x=zeros（n,1）;
fori=1:
nx1（i）=i;
endA=zeros（n,n）;fori=1:
1:
n
forj=1:
1:
nifi>jA（i,j）=-1;
elseif （ielseA（i,j）=1;
endend
endb=A*x1;
[P,L,U]=PLU（A）;
b1=b;b=P*b1;
y=zeros（n,1）;
y
（1）=b
（1）;
fori=2:
n
y（i）=b（i）-sum（L（i,1:
i-1）'.*y（1:
i-1））;end
%yx（n）=y（n）/U（n,n）;fori=n-1:
-1:
1
x（i）=（y（i）-sum（U（i,i+1:
n）'.*x（i+1:
n）））/U（i,i）;end
%x
d=norm（x1-x）;%计算算法误差end
运行：
cleard=zeros（26,1）;forn=5:
30
d（n-4）=PLU1（n）;
endd1=d';d1
d1=
Columns1through5
0 0 0 0 0
Columns6through10
0 0 0 0 0
Columns11through15
0 0 0 0 0
Columns16through20
0 0 0 0 0
Columns21through25
0 0 0 0 0
Column26
0
由程序计算结果可知，PLU分解计算程序计算误差非常小，基本无误差计算
（2）
%由PLU分解求逆矩阵程序
functionB=PLU2（n）x=zeros（n,1）;
A=zeros（n,n）;
B=zeros（n,n）;fori=1:
1:
n
forj=1:
1:
nifi>jA（i,j）=-1;
elseif （iA（i,j）=0;
elseA（i,j）=1;
endend
ende=ones（n,n）;fork=1:
n
b（k）=e（:
k）;
[P,L,U]=PLU（A）;
b1=b（k）;b=P*b1;y=zeros（n,1）;
y
（1）=b
（1）;
fori=2:
n
y（i）=b（i）-sum（L（i,1:
i-1）'.*y（1:
i-1））;end
%yx（n）=y（n）/U（n,n）;fori=n-1:
-1:
1
x（i）=（y（i）-sum（U（i,i+1:
n）'.*x（i+1:
n）））/U（i,i）;end
B（:
k）=x;
end
P=A*B%验证是否正确
%x
%d=norm（x1-x）;%计算算法误差end
clear
forn=5:
30
B{n-4}=zeros（n）;
B{n-4}=PLU2（n）
end
求出的n从5到30的逆矩阵存储在B中：
B=
Columns1through5
[5x5double] [6x6double] [7x7double] [8x8double] [9x9
double]
Columns6through9
[10x10double] [11x11double] [12x12double] [13x13double]
Columns10through13
[14x14double] [15x15double] [16x16double] [17x17double]
Columns14through17
[18x18double] [19x19double] [20x20double] [21x21double]
Columns18through21
[22x22double] [23x23double] [24x24double] [25x25double]
Columns22through25
[26x26double] [27x27double] [28x28double] [29x29double]
Column26
[30x30double]
5.编制计算对称正定阵的Cholesky分解的通用程序，并用你编制的程序计算
Ax=b，其中A=A（i,j）ÎRn´n,a
= 1
ij
i+j-1
。
b可以由你自己取定，对n从10
到20验证程序的可靠性。
%Cholesky分解的通用程序
functionL=Cholesk（A）
%如果是正定矩阵，可以进行下面分解操作L（1,1）=sqrt（A（1,1））;%确定第一列元素n=length（A）;
fori=2:
n
L（i,1）=A（i,1）/L（1,1）;
endsum1=0;
forj=2:
n-1%确定第j列元素
fork=1:
j-1
sum1=sum1+L（j,k）^2;
end
L（j,j）=sqrt（A（j,j）-sum1）;sum1=0;
fori=j+1:
n
fork=1:
j-1
sum1=sum1+L（i,k）*L（j,k）;
end
L（i,j）=（A（i,j）-sum1）/L（j,j）;
endsum1=0;
end
fork=1:
n-1%确定第n列元素
sum1=sum1+L（n,k）^2;
end
L（n,n）=sqrt（A（n,n）-sum1）;
P=L*L';%验证分解是否正确end
%Cholesky分解的通用程序用于生成A,b矩阵的函数程序
function[x1,b,A]=Ch（n）x1=zeros（n,1）;
fork=1:
n
x1（k）=1;
endA=zeros（n,n）;fori=1:
n
forj=1:
n
A（i,j）=1/（i+j-1）;
endendb=A*x1;end
%计算Cholesky分解的通用程序的误差程序
functione=Cho（n）
[x1,b,A]=Ch（n）;
[m,n]=size（A）;
ifm~=n%判断输入的矩阵是不是方阵e=inf;
return;end
d=eig（A）;%根据方阵的特征值判定是不是正定矩阵
fori=1:
n
ifd（i）<=0e=inf;
return;else
break;
end
end
%如果是正定矩阵，可以进行下面分解操作L=Cholesk（A）;
y=zeros（n,1）;y
（1）=b
（1）/L（1,1）;
fori=2:
n
y（i）=（b（i）-sum（L（i,1:
i-1）'.*y（1:
i-1）））/L（i,i）;end
x（n）=y（n）/L（n,n）;fori=n-1:
-1:
1
x（i）=（y（i）-sum（L（i+1:
n,i）'.*x（i+1:
n）））/L（i,i）;
endx2=x';x=x2;
%A
%x1
%x
e=norm（x1-x）;%计算算法误差end
运行：
>>cleard=zeros（19,1）;forn=2:
20
d（n-1）=Cho（n）;end
d
d=
1.0e+17*
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000042
0.000000000000000
0.000000001856285
0.000000000000000
0.000046586561891
0.000000000007943
1.394903786960288
0.000000118884207
Inf0.002275657682969
InfInfInfInfInfInf
由上面程序运算结果是n从2~20的的程序误差结果可知，低阶是程序可靠性较好，随着n的增大程序的误差总体越来越大，（运行结果后面输出为inf（程序设定）的n阶矩阵不是正定矩阵）程序的可靠性变差。
6.
（1）编制程序House（x），其作用是对输入的向量x，输出单位向量u使得
x
（I-2uuT）x=
2e1。
（2）编制Householder变换阵H=I-2uuTÎRn´n乘以AÎRn´m的程序HA，注意，你的程序并不显式的计算出H。
（3）考虑矩阵
3
2
4
3
æ 1 2 ö
ç-1 3 ÷
A=ç ÷
ç-10 2 -3 7÷
ç
è 0 2
÷
7 5/2ø
用你编制的程序计算H使得HA的第一列为ae1的形式，并将HA的结果显示。
（1）
functionu=house（x）e=eye（length（x））;w=x-norm（x）*e（:
1）;u=w./norm（w）;
升级会员 