级矩阵与数值分析上机作业.docx
《级矩阵与数值分析上机作业.docx》由会员分享，可在线阅读，更多相关《级矩阵与数值分析上机作业.docx（53页珍藏版）》请在冰豆网上搜索。
级矩阵与数值分析上机作业
2016级矩阵与数值分析上机作业
学生班级：
学生姓名：
任课教师：
所在学院：
电子信息与电气学部
学生学号：
使用软件：
MATLAB
2016年12月10号
1.考虑计算给定向量的范数：
输入向量
，输出
，
，
请编制一个通用程序，并用你编制的程序计算如下向量的范数：
，
对n=10，100，1000甚至更大的n计算其范数，你会发现什么结果？
你能否修改你的程序使得计算结果相对精确呢？
（1）计算范数的程序
function[Nor1,Nor2,Nor3]=Nor（a）
b=length（a）;
formatlong
Nor11=0;
formatlong
Nor21=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））;
end
i=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;
end
end
functionb=bfun（n）
b=zeros（n）;
b
（1）=1;
fori=1:
1:
n;
b（i）=i;
end
end
运行：
>>n=10;
>>a=afun（n）;
>>b=bfun（n）;
>>[f1,f2,f3]=Nor（a）;
>>[h1,h2,h3]=Nor（b）;
f1=2.928968253968254f2=1.244896674895769f3=1
h1=55h2=19.621416870348583h3=10
>>n=100;
a=afun（n）;
b=bfun（n）;
[f1,f2,f3]=Nor（a）;
>>[h1,h2,h3]=Nor（b）;
f1=5.9621f2=1.278664889713052f3=1
h1=5050h2=5.816786054171153e+02h3=100
>>n=1000;
a=afun（n）;
b=bfun（n）;
[f1,f2,f3]=Nor（a）;
>>[h1,h2,h3]=Nor（b）;
f1=7.485470860550343f2=1.282160117411847f3=1
h1=500500h2=1.827111107732642e+04h3=1000
上述结果由于MATLAB的高精度运算较为准确，但是由于当n非常大时可能会出现了大数吃小数的现象使误差增大，而y向量的范数结果较为准确是由于y向量是按从小到大排列。
改进范数计算程序使输入序列按从小到大排列再计算其范数：
function[Nor1,Nor2,Nor3]=Nor（c）
a=sort（c）;b=length（a）;
formatlong
Nor11=0;
formatlong
Nor21=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））;
end
i=i+1;
end
Nor1=Nor11;%计算1范数
Nor2=sqrt（Nor21）;%计算2范数
Nor3=t;%计算无穷范数
end
2.考虑
，其中定义
，此时
是连续函数。
用此公式计算当
时的函数值，画出图像。
另一方面，考虑下面算法：
用此算法计算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
运行：
clc
clear
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.首先编制一个利用秦九韶算法计算一个多项式在给定点的函数值的通用程序，你的程序包括输入多项式的系数以及给定点，输出函数值。
利用你编制的程序计算
在x=2邻域附近的值。
画出p（x）在x∈[1.95,20.5]上的图像。
秦九韶算法程序：
functionp=p（a,x）
n=length（a）;
p=a（n）;
fori=n:
-1:
2
pi=p.*x+a（i-1）;
p=pi;
end
end
运行：
clc
clear
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）考虑：
自己取定
，并计算b=Ax，然后利用你编制的不选主元和列主元的Gausss消去法求解该方程组，记你算的解为
，对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误差分析程序
functiond=zlu1（n）
x1=zeros（n,1）;
x=zeros（n,1）;
fori=1:
n
x1（i）=i;
end
A=zeros（n,n）;
fori=1:
1:
n
forj=1:
1:
n
ifi>j
A（i,j）=-1;
elseif（iA（i,j）=0;
else
A（i,j）=1;
end
end
end
b=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
%y
x（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
运行：
clear
d=zeros（26,1）;
forn=5:
30
d（n-4）=zlu1（n）;
end
d1=d';
d1
d1=
Columns1through9
000000000
Columns10through18
000000000
Columns19through26
00000000
由程序计算结果可知，LU分解计算程序计算误差非常小，基本无误差计算
2.PLU分解
%PLU分解,PA=LU
function[P,L,U]=PLU（A）
formatrat
B=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）==0
fork=i+1:
m
ifB（k,i）~=0
x=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）;
end
B（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:
n
x1（i）=i;
end
A=zeros（n,n）;
fori=1:
1:
n
forj=1:
1:
n
ifi>j
A（i,j）=-1;
elseif（iA（i,j）=0;
else
A（i,j）=1;
end
end
end
b=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
%y
x（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
运行：
clear
d=zeros（26,1）;
forn=5:
30
d（n-4）=PLU1（n）;
end
d1=d';
d1
d1=
Columns1through5
00000
Columns6through10
00000
Columns11through15
00000
Columns16through20
00000
Columns21through25
00000
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:
n
ifi>j
A（i,j）=-1;
elseif（iA（i,j）=0;
else
A（i,j）=1;
end
end
end
e=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
%y
x（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][9x9double]
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，其中
。
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）;
end
sum1=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）;
end
sum1=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;
end
A=zeros（n,n）;
fori=1:
n
forj=1:
n
A（i,j）=1/（i+j-1）;
end
end
b=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）<=0
e=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）;
end
x2=x';
x=x2;
%A
%x1
%x
e=norm（x1-x）;%计算算法误差
end
运行：
>>clear
d=zeros（19,1）;
forn=2:
20
d（n-1）=Cho（n）;
end
d
d=
1.0e+17*
0.0000
0.0000
0.0000
0.0042
0.0000
0.6285
0.0000
0.1891
0.7943
1.3949
0.4207
Inf
0.2969
Inf
Inf
Inf
Inf
Inf
Inf
由上面程序运算结果是n从2~20的的程序误差结果可知，低阶是程序可靠性较好，随着n的增大程序的误差总体越来越大，（运行结果后面输出为inf（程序设定）的n阶矩阵不是正定矩阵）程序的可靠性变差。
6.
（1）编制程序
，其作用是对输入的向量x，输出单位向量u使得
。
（2）编制Householder变换阵
乘以
的程序HA，注意，你的程序并不显式的计算出
。
（3）考虑矩阵
用你编制的程序计算
使得HA的第一列为
的形式，并将HA的结果显示。
（1）
functionu=house（x）
e=eye（length（x））;
w=x-norm（x）*e（:
1）;
u=w./norm（w）;
end
（2）
function[H,HA]=HA（x,A）
u=house（x）;
H=eye（length（u））-2*u*u';
HA=H*A;
end
（3）
clear
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,2.5];
x=A（:
1）;
[H,HA]=HA（x,A）
H=
0.2500-0.2500-0.5000-0.79060
-0.25000.9167-0.1667-0.26350
-0.5000-0.16670.6667-0.52700
-0.7906-0.2635-0.52700.16670
00001.0000
HA=
4.0000-2.83111.4090-6.5378
-0.00001.38960.8839-1.7805
-0.0000-1.22081.6576-3.8836
-0.0000-3.0925-4.6770-4.1078
02.00007.00002.5000
7.用Jacobi和Gauss-Seidel迭代求解下面的方程组，输出迭代每一步的误差
：
假设迭代法计算停止的条件为：
（1）Jacobi迭代法
clc;
clear;
A=[5,-1,-3;-1,2,4;-3,4,15];
b=[-2,1,10];
delta=10^（-4）;%误差
n=length（A）;
k=0;
x=zeros（n,1）;
y=zeros（n,1）;
while1
fori=1:
n
y（i）=b（i）;
forj=1:
n
ifj~=i
y（i）=y（i）-A（i,j）*x（j）;
end
end
y（i）=y（i）/A（i,i）;%迭代出的x
end
d=norm（y-x）;
ifdbreak;
end
x=y;
k=k+1;
B（1,k）=k;
B（2,k）=d;
end
B
y
运行结果：
（B的第一行是迭代次数k，第二行是每一步迭代误差
）
B=
Columns1through8
1.00002.00003.00004.00005.00006.00007.00008.0000
0.92441.62680.95171.33830.95421.14900.92611.0146
Columns9through16
9.000010.000011.000012.000013.000014.000015.000016.0000
0.87940.91140.82420.82690.76650.75470.70970.6912
Columns17through24
17.000018.000019.000020
升级会员 