大连理工大学矩阵分析matlab上机作业

矩阵与数值分析上机实验报告学生姓名：高明星班级：矩阵与数值分析4班学号：31603026专业：船舶与海洋工程授课老师：董波大连理工大学DalianUniversityofTechnology1计算给定向量的范数题目：考虑计算给定向量的范数：输入向量x=(𝑥1,𝑥2,…,𝑥𝑛)𝑇,输出‖𝑥‖1,‖𝑥‖2,‖𝑥‖∞。请编制一个通用程序，并用你编制的程序计算如下向量的范数：x=(1,12,13,…,1𝑛)𝑇,𝑦=(1,2,…,𝑛)𝑇.对n=10，100，1000甚至更大的n计算其范数，你会发现什么结果?你能否修改你的程序使得计算结果相对精确呢？1.1源代码function[z1,z2]=norm_vector(n)%向量z1的值为向量x的是三种范数，向量z2的值为向量y的三种范数，n为输入参数x=zeros(n,1);%为列向量x预分配存储空间y=1:n;%定义行向量yy=y';%把行向量y改成列向量fori=1:nx(i)=1/i;%按要求给向量x赋值，其值递减endnormx1=norm(x,1);%求解向量x的1范数normx1normx2=norm(x,2);%求解向量x的2范数normx2normxinf=norm(x,inf);%求解向量x的无穷范数normxinfnormy1=norm(y,1);%求解向量y的1范数normy1normy2=norm(y,2);%求解向量y的2范数normy2normyinf=norm(y,inf);%求解向量y的无穷范数normyinfz1=[normx1,normx2,normxinf];z2=[normy1,normy2,normyinf];end1.2函数调用及运行结果n=10,[z1z2]=norm_vector(10)n=100,[z1z2]=norm_vector(100)n=1000,[z1z2]=norm_vector(1000)n=10z1=2.9289682539682541.2448966748957691.000000000000000z2=55.00000000000000019.62141687034858310.000000000000000n=100z1=5.1873775176396211.2786648897130521.000000000000000z2=1.0e+03*5.0500000000000000.5816786054171150.100000000000000n=1000z1=7.4854708605503441.2821601174118471.000000000000000z2=1.0e+05*5.0050000000000000.1827111107732640.010000000000000将以上数据汇总为表1-1表1-1向量x,y各范数汇总n=10n=100n=10001x2.9289682539682545.1873775176396217.4854708605503442x1.2448966748957691.2786648897130521.282160117411847x1.0000000000000001.0000000000000001.0000000000000001y55.0000000000000005050.0000000000005.00500.00000000002y19.621416870348583581.67860541711518271.1110773264y10.0000000000000001000.000000000001000.0000000000结果分析：该源代码在求解向量x的1范数和2范数时存在“大数吃小数”的问题，会降低其精度，带来误差，当n越大时，误差会越大，因此可逆序求和。修改后的程序代码为:functionz1=norm_vector2(n)%向量z1的值为向量x的是三种范数，向量z2的值为向量y的三种范数%n为输入参数x=zeros(n,1);%为列向量x预分配存储空间fori=1:nx(n-i+1)=1/(i);%按升序给向量x赋值，其值递增endnormx1=norm(x,1);%求解向量x的1范数normx1normx2=norm(x,2);%求解向量x的1范数normx2normxinf=norm(x,inf);%求解向量x的1范数normxinfz1=[normx1,normx2,normxinf];end取n=106、108分别调用上述两种程序，所得结果如表1-2，通过比较可以发现，两种方法所得向量x的1范数和2范数相差不大，但随着n的增大，两者之间的差值在变大。表1-2修改前后向量x的1范数、2范数对比n=106n=108修改前修改后修改前修改后1x14.39272672286574414.39272672286572318.99789641385351818.9978964138539122x1.2825494403136161.2825494403135991.2825498266479071.2825498262633802计算函数的值题目：考虑y=f(x)=ln⁡(1+𝑥)𝑥,其中定义f(0)=1，此时f(x)是连续函数。用此公式计算当x∈[−10−15,10−15]时的函数值，画出图像。另一方面，考虑下面算法：d=1+xifd=1theny=1elsey=lnd/(d-1)endif用此算法计算x∈[−10−15,10−15]时的函数值，画出图像。比较一下发生了什么？2.1源程序及运行结果clcx=(-1:0.01:1)*1E-15;y1=log(1+x)./x;d=1+x;fori=1:length(x)ifd(i)==1y2(i)=1;elsey2(i)=log(d(i))/(d(i)-1);endendplot(x,y1,'--r',x,y2,'m');xlabel('x');ylabel('y');legend('修改前','修改后');图2-1修改前后y=ln(1+x)/x函数值对比结果分析：图2-1为采用两种算法计算函数y=ln(1+x)/x在区间[10-15,1015]上的图像对比，可以看出第一种算法所画图像在y=1处上下震荡，呈现出数值的不稳定，这主要是由于小数做除数损失了有效数字，而第二种算法避免了有效数字的损失，所画出的图像在x=0附近均为1，数值很稳定。3秦九韶算法题目：首先编制一个利用秦九韶算法计算一个多项式在给定点的函数值的通用程序，你的程序包括输入多项式的系数以及给定点，输出函数值。利用你编制的程序计算p(x)=(𝑥−2)9=𝑥9−18𝑥8+144𝑥7−672𝑥6+2016𝑥5−4032𝑥4+5376𝑥3−4608𝑥2−512在x=2邻域附近的值。画出p(x)在x∈[1.95,20.5]上的图像。-1-0.500.51x10-1500.511.52xy修改前修改后3.1源程序用秦九韶算法计算一个多项式在给定点的函数值的通用程序为：functiony=QJS(a,x)%QJS()为用秦九韶算法计算多项式值的函数%a为多项式的系数向量，x为自变量，y为因变量n=length(a);y=a(1);fori=2:ny=y*x+a(i);endend3.2函数调用及运行结果函数调用：clca=[1-18144-6722016-40325376-46082304-512];x=2;y=QJS(a,x)z=linspace(1.95,2.05,100);y1=zeros(1,length(z));y2=zeros(1,length(z));fori=1:length(z)y1(i)=QJS(a,z(i));y2(i)=(z(i)-2)^9;i=i+1;endplot(z,y1);holdon;plot(z,y2,'or');xlabel('x');ylabel('y');legend('秦九韶公式','y=(x-2)^9');运行结果：y=0图3-1y=(x-2)9的图像结果分析：当x=2时，通过调用秦九韶函数QJS()计算多项式y=(x-2)9，得到y=0，初步验证了该通用程序的正确性。图3-1为多项式y=(x-2)9在区间[1.95,20.5]上用该程序求得的解与真实解绘制的对比图，可以发现，在零点领域内，用秦九韶求得的解是不稳定的，这主要是由于相近数相减损失了有效数字。4LU分解和PLU分解题目：编制计算给定矩阵A的LU分解和PLU分解的通用程序，然后用你编制的程序完成下面两个计算任务：(1)考虑A=[10⋯−1⋱⁡⁡⋱⋮⋱⋱⁡⁡⁡⁡⁡⁡0⋮0⁡⁡⁡⁡⁡⁡⁡1⋮1−1⋯−1⁡⁡⁡⁡⁡⁡⁡11−1⋯−1⁡⁡⁡⁡−11]∈𝑅𝑛×𝑛,自己取定x∈𝑅𝑛，并计算b=Ax。然后用你编制的不选主元和列主元的Gauss消去法求解该方程组，记你计算出的解为𝑥̂。对n从5到30估计计算解的精度。(2)对n从5到30计算其逆矩阵。4.1源程序生成给定矩阵A的函数matrix(n)定义：functionA=matrix(n)%n为矩阵阶数，A为按要求生成的给定矩阵L=-tril(ones(n),-1);D=diag(ones(n,1));U=zeros(n);1.9522.05-1.5-1-0.500.511.5x10-11xy秦九韶公式y=(x-2)9fori=1:n-1U(i,n)=1;endA=L+D+U;endLU分解函数LU(A)定义：function[L,U]=LU(A)%求矩阵A的LU分解，返回单位下三角矩阵L，上三角矩阵Un=length(A);L=eye(n);U=zeros(n);U(1,:)=A(1,:);L(:,1)=A(:,1)/U(1,1);fori=2:nforj=i:nU(i,j)=A(i,j)-L(i,1:i-1)*U(1:i-1,j);%Doolittle分解计算上三角矩阵的公式L(j,i)=(A(j,i)-L(j,1:i-1)*U(1:i-1,i))/U(i,i);%Doolittle分解计算下三角矩阵的公式endendendPLU分解函数PLU(A)定义：function[P,L,U]=PLU(A)%求矩阵A的LU分解，返回单位下三角矩阵L，上三角矩阵U，置换矩阵Pn=length(A);P=eye(n);L=eye(n);index=1;%用于显示列主元所在行号forj=1:n-1[mainelement,index]=max(abs(A(j:n,j)));index=index+j-1;%将mainelement在A(j:n,j)中所在行号转化为在A中的行号ifindexj%列主元所在行不是当前行，需要行交换temp1=A(j,:);A(j,:)=A(index,:);A(index,:)=temp1;%矩阵A行交换temp2=P(j,:);P(j,:)=P(index,:);P(index,:)=temp2;%置换矩阵P行交换endfori=j+1:nL(i,j)=A(i,j)/A(j,j);%求下三角矩阵i列元素A(i,j:n)=A(i,j:n)-L(i,j)*A(j,j:n);endendU=A;endLU分解求解Ax=b函数LUsolve(A,b)的定义：functionx=LUsolve(A,b)%用LU分解求解Ax=b[L,U]=LU(A);n=length(A);x=zeros(n,1);y=zeros(n,1);y(1)=b(1)/L(1,1);fori=2:ny(i)=(b(i)-sum(L(i,1:i-1)*y(1:i-1)))/L(i,i);%求解Ly=bendx(n)=y(n)/U(n,n);fori=n-1:-1:1x(i)=(y(i)-sum(U(i,i+1:n)*x(i+1:n)))/U(i,i);%求解Ux=yendendPLU分解求解Ax=b函数PLUsolve(A,b)的定义：funct