级矩阵与数值分析上机作业.docx
- 文档编号:18102753
- 上传时间:2023-08-13
- 格式:DOCX
- 页数:61
- 大小:830.49KB
级矩阵与数值分析上机作业.docx
《级矩阵与数值分析上机作业.docx》由会员分享,可在线阅读,更多相关《级矩阵与数值分析上机作业.docx(61页珍藏版)》请在冰点文库上搜索。
级矩阵与数值分析上机作业
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.187377517639621f2=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(i A(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(i A(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(i A(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.000000000000000 0.000000000000000 0.000000000000000 0.000000000000042 0.000000000000000 0.000000001856285 0.000000000000000 0.000046586561891 0.000000000007943 1.394903786960288 0.000000118884207 Inf 0.002275657682969 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); ifd break; 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.0000
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 矩阵 数值 分析 上机 作业