数值线性代数第二版徐树方高立张平文上机习题第二章实验报告

合集下载

线性代数第二章习题部分答案

线性代数第二章习题部分答案

线性代数第二章习题部分答案第二章向量组的线性相关性§2-1 §2-2 维向量,线性相关与线性无关(一)一、填空题1. 设3 α1α +2 α2+α =5 α3+α , 其中α1=(2,5,1,3)T,α2=(10,1,5,10)T, α3=(4,1,1,1)T, 则α= (1,2,3,4)T .2. 设α1=(1,1,1)T, α2=(2,1,1)T,α3=(0,2,4)T,则线性组合α13α2+α3= (5,0,2)T .3. 设矩阵A= 5 ,设βi为矩阵A的第i个列向量,则2β1+β2β3= (2,8,2)T .二、试确定下列向量组的线性相关性1. α1=(2,1,0)T, α2=(1,2,1)T, α3=(1,1,1)T解:设k1α1+k2α2+k3α3=0,则k1 210 +k2 121 +k3 111 = 000即2k1+k2+k3=0k1+2k2+k3=0k2+k3=0 k1+2k2+k3=03k2k3=0k2+k3=0 k1+2k2+k3=0k2+k3=0k3=0 k1=k2=k3=0,线性无关。

2. α1=(1,1,2)T, α2=(0,0,0)T, α3=(1,4,3)T线性相关三、设有向量组α1=(1,1,0)T, α2=(1,3,1)T, α3=(5,3,t)T,问t取何值时该向量组线性相关。

解:设k1α1+k2α2+k3α3=0,则k1 110 +k2 131 +k3 53t =0即k1+k2+5k3=0k1+3k23k3=0k2+tk3=0 k1+k2+5k3=0k24k3=0k2+tk3=0k1+k2+5k3=0k1+3k23k3=0(t4)k3=0所以,t=4, 线性相关; t≠4, 线性无关四、设a1,a2线性无关,a1+b,a2+b线性相关,求向量b用a1,a2线性表示的表示式。

解:因为a1+b,a2+b线性相关,所以存在不全为零的k1,k2,使得k1(a1+b)+k2(a2+b)=0, 即(k1+k2)b=k1a1k2a2.又因为a1,a2线性无关,所以k1+k2≠0,于是,b=k1k1+k2a1k2k1+k2a2.五、已知向量组α1,α2,,α2n,令β1=α1+α2,β2=α2+α3,,β2n=α2n+α1,求证向量组β1,β2,,β2n线性相关。

数值代数上机实验报告

数值代数上机实验报告

数值代数上机实验报告试验项目名称:平方根法与改进平方根法实验内容:先用你熟悉的计算机语言将平方根法和改进平方根法编写成通用的子程序,然后用你编写的程序求解对称正定方程组Ax=b,其中,A=[101 10 1…1 10 11 10]100*100b随机生成,比较计算结果,评论方法优劣。

实验要求:平方根法与改进的平方根的解法步骤;存储单元,变量名称说明;系数矩阵与右端项的生成;结果分析。

实验报告姓名:罗胜利班级:信息与计算科学0802 学号:u200810087 实验一、平方根法与改进平方根法先用你所熟悉的计算机语言将平方根法和改进的平方根法编成通用的子程序,然后用你编写的程序求解对称正定方程组AX=b,其中系数矩阵为40阶Hilbert矩阵,即系数矩阵A的第i行第j列元素为=,向量b的第i个分量为=.平方根法函数程序如下:function [x,b]=pingfanggenfa(A,b)n=size(A);n=n(1);x=A^-1*b; %矩阵求解disp('Matlab自带解即为x');for k=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endend %Cholesky分解for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n); %前代法A=A';for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1); %回代法disp('平方根法的解即为b');endfunction [x]=ave(A,b,n) %用改进平方根法求解Ax=b L=zeros(n,n); %L为n*n矩阵D=diag(n,0); %D为n*n的主对角矩阵S=L*D;for i=1:n %L的主对角元素均为1L(i,i)=1;for i=1:nfor j=1:n %验证A是否为对称正定矩阵if (eig(A)<=0)|(A(i,j)~=A(j,i)) %A的特征值小于0或A非对称时,输出wrong disp('wrong');break;endendendD(1,1)=A(1,1); %将A分解使得A=LDL Tfor i=2:nfor j=1:i-1S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);endD(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');endy=zeros(n,1); % x,y为n*1阶矩阵x=zeros(n,1);for i=1:ny(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i); %通过LDy=b解得y的值endfor i=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n)); %通过L T x=y解得x的值改进平方根法函数程序如下:function b=gaijinpinfanggenfa(A,b)n=size(A);n=n(1);v=zeros(n,1);for j=1:nfor i=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);end %LDL'分解B=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;EndA=tril(A); %得到L和Dfor j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n); %前代法A=D*(A');for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1); %回代法disp('改进平方根法解得的解即为b');end调用函数解题:clear;clc;n=input('请输入矩阵维数:');b=zeros(n,1);A=zeros(n);for i=1:nfor j=1:nA(i,j)=1/(i+j-1);b(i)=b(i)+1/(i+j-1);endend %生成hilbert矩阵[x,b]=pingfanggenfa(A,b) b=gaijinpinfanggenfa(A,b)运行结果:请输入矩阵维数:40Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.570692e-020. > In pingfanggenfa at 4In qiujie at 10Matlab自带解即为x平方根法的解即为bx =1.60358.96850.85621.01950.9375-50.2500-3.0000-16.000024.0000-49.5000-30.000039.000022.0000-64.0000-12.00002.000010.2500-10.5000-1.0000-10.875083.000046.0000-98.0000-69.000068.000021.0000-50.7188-8.7500-8.0000 112.0000 6.0000 -68.7500 22.000044.0000 -28.0000 8.0000 -44.000012.0000b =1.0e+007 *0.0000-0.00000.0001-0.0004-0.00140.0424-0.29801.1419-2.73354.2539-4.30182.7733-1.19890.5406-0.36880.32850.4621-0.25130.05650.0000-0.00510.0071-0.0027-0.0031-0.00190.00090.0002-0.0002-0.00060.00040.0001-0.00020.00010.0000-0.00000.0000-0.0000-0.0000改进平方根法解得的解即为bb =1.0e+024 *0.0000-0.00000.0001-0.0012-0.0954 0.4208 -1.2101 2.0624 -1.0394 -3.3343 6.2567 -0.2463 -7.45942.80303.6990 0.7277 -1.7484 -0.4854 -3.6010 0.2532 5.1862 1.4410 0.8738 -4.5654 1.0422 4.0920 -2.7764 -2.2148 -0.8953 0.3665 4.8967 1.0416 0.1281-1.1902-2.83348.4610-3.6008实验二、利用QR分解解线性方程组:利用QR分解解线性方程组Ax=b,其中A=[16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12];b=[32 26 38 30];求解程序如下:定义house函数:function [v,B]=house(x)n=length(x);y=norm(x,inf);x=x/y;Q=x(2:n)'*x(2:n);v(1)=1;v(2:n)=x(2:n);if n==1B=0;elsea=sqrt(x(1)^2+Q);if x(1)<=0v(1)=x(1)-a;elsev(1)=-Q/(x(1)+a);endB=2*v(1)^2/(Q+v(1)^2);endend进行QR分解:clear;clc;A=[16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12]; b=[32 26 38 30];b=b';x=size(A);m=x(1);n=x(2);d=zeros(n,1);for j=1:n[v,B]=house(A(j:m,j));A(j:m,j:n)=(eye(m-j+1)-B*(v')*v)*A(j:m,j:n); d(j)=B;if j<m< p="">A(j+1:m,j)=v(2:m-j+1);endend %QR分解R=triu(A); %得到R D=A;I=eye(m,n);Q=I;for i=1:nD(i,i)=1;endH=tril(D);M=H';for i=1:nN=I-d(i)*H(1:m,i)*M(i,1:m);Q=Q*N;end %得到Qb=(Q')*b; %Q是正交阵for j=n:-1:2b(j)=b(j)/R(j,j);b(1:j-1)=b(1:j-1)-b(j)*R(1:j-1,j);endb(1)=b(1)/R(1,1); %回带法运行结果如下:R =18.7617 9.8072 15.7769 11.08640 9.9909 9.3358 7.53410 0 5.9945 9.80130 0 0 -0.5126Q =0.8528 -0.4368 -0.2297 -0.17090.2132 0.7916 -0.4594 -0.34170.4264 0.3822 0.2844 0.76890.2132 0.1911 0.8095 -0.5126b=1.000000000000001.000000000000010.9999999999999881.00000000000001实验三、Newton下山法解非线性方程组:3x-cos(yz)-=0,-81+sinz+1.06=0,exp(-xy)+20z+=0;要求满足数值解=满足或.定义所求方程组的函数:Newtonfun.mfunction F = Newtonfun(X)F(1,1)=3*X(1)-cos(X(2)*X(3))-1/2;F(2,1)=X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;F(3,1)=exp(-X(1)*X(2))+20*X(3)+(10*pi-3)/3;End向量求导:Xiangliangqiudao.mfunction J=xiangliangqiudao()syms x y zX=[x,y,z];F=[3*X(1)-cos(X(2)*X(3))-1/2;X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;exp(-X(1)*X(2))+20*X(3)+(10*pi-3)/3];J=jacobian(F,[x y z]);End代值函数:Jacobi.mfunction F=Jacobi(x)F=[ 3,x(3)*sin(x(2)*x(3)), x(2)*sin(x(2)*x(3));2*x(1), -162*x(2)-81/5,cos(x(3));-x(2)/exp(x(1)*x(2)),-x(1)/exp(x(1)*x(2)),20];End方程组求解:format long; %数据表示为双精度型X1=[0,0,0]';eps=10^(-8);k=1;i=1;X2=X1-Jacobi(X1)^(-1)*Newtonfun(X1);while (norm(subs(X2-X1,pi,3.1415926),2)>=eps)&&(norm(Newtonfun(X1),2)>=eps) if norm(Newtonfun(X2),2)<="" p="">X1=X2;B=inv(Jacobi(X2));C=Newtonfun(X2);X2=X2-B*C;i=i+1;elsev=1/(2^k); %引入下山因子X1=X2;B=inv(Jacobi(X2));C=Newtonfun(X2);X2=X2-v*B*C;k=k+1;endendj=i+k-1 %迭代次数X=X2 %输出结果运行结果如下:j =5X =0.500000000000000 -0.000000000000000 -0.523598775598299</m<>。

数值线性代数第二版徐树方高立张平文上机习题第一章实验报告(供参考)

数值线性代数第二版徐树方高立张平文上机习题第一章实验报告(供参考)

上机习题1.先用你所熟悉的的计算机语言将不选主元和列主元Gauss 消去法编写成通用的子程序;然后用你编写的程序求解84阶方程组;最后将你的计算结果与方程的精确解进行比较,并就此谈谈你对Gauss 消去法的看法。

Sol :(1)先用matlab 将不选主元和列主元Gauss 消去法编写成通用的子程序,得到P U L ,,: 不选主元Gauss 消去法:[])(,A GaussLA U L =得到U L ,满足LU A =列主元Gauss 消去法:[])(,,A GaussCol P U L =得到P U L ,,满足LU PA =(2)用前代法解()Pb or b Ly =,得y用回代法解y Ux =,得x求解程序为()P U L b A Gauss x ,,,,=(P 可缺省,缺省时默认为单位矩阵)(3)计算脚本为ex1_1代码%算法(计算三角分解:Gauss 消去法)function [L,U]=GaussLA(A)n=length(A);for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endU=triu(A);L=tril(A);L=L-diag(diag(L))+diag(ones(1,n));end%算法计算列主元三角分解:列主元Gauss消去法)function [L,U,P]=GaussCol(A)n=length(A);for k=1:n-1[s,t]=max(abs(A(k:n,k)));p=t+k-1;temp=A(k,1:n);A(k,1:n)=A(p,1:n);A(p,1:n)=temp;u(k)=p;if A(k,k)~=0A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); elsebreak;endendL=tril(A);U=triu(A);L=L-diag(diag(L))+diag(ones(1,n));P=eye(n);for i=1:n-1temp=P(i,:);P(i,:)=P(u(i),:);P(u(i),:)=temp;endend%高斯消去法解线性方程组function x=Gauss(A,b,L,U,P)if nargin<5P=eye(length(A));endn=length(A);b=P*b;for j=1:n-1b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j); endb(n)=b(n)/L(n,n);y=b;for j=n:-1:2y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endy(1)=y(1)/U(1,1);x=y;endex1_1clc;clear;%第一题A=6*eye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1);b=[7;15*ones(82,1);14];%不选主元Gauss消去法[L,U]=GaussLA(A);x1_1=Gauss(A,b,L,U);%列主元Gauss消去法[L,U,P]=GaussCol(A);x1_2=Gauss(A,b,L,U,P);%解的比较subplot(1,3,1);plot(1:84,x1_1,'o-');title('Gauss');subplot(1,3,2);plot(1:84,x1_2,'.-');title('PGauss');subplot(1,3,3);plot(1:84,ones(1,84),'*-');title('精确解');结果为(其中Gauss表示不选主元的Gauss消去法,PGauss表示列主元Gauss 消去法,精确解为[]'⨯8411,,1 ):-6-4-202468Gauss050100PGauss 00.20.40.60.811.21.41.61.82精确解由图,显然列主元消去法与精确解更为接近,不选主元的Gauss 消去法误差比列主元消去法大,且不如列主元消去法稳定。

2021年数值线性代数第二版徐树方高立张平文上机习题实验报告2

2021年数值线性代数第二版徐树方高立张平文上机习题实验报告2

第四章上机习题1考虑两点边值问题⎪⎩⎪⎨⎧==<<=+.1)1(,0)0(10 ,22y y a a dx dy dx y d ε 轻易知道它正确解为ax e e ay x +---=--)1(111εε为了把微分方程离散化, 把[0,1]区间n 等分, 令h=1/n,1,,1,-==n i ih x i得到差分方程,21211a hy y h y y y i i i i i =-++-++-ε简化为 ,)2()(211ah y y h y h i i i =++-+-+εεε从而离散化后得到线性方程组系数矩阵为⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡+-++-++-++-=)2()2()2()2(h h h h h h h A εεεεεεεεεε 对,100,2/1,1===n a ε分别用Jacobi 迭代法, G-S 迭代法和SOR 迭代法求线性方程组解, 要求有4位有效数字, 然后比较与正确解得误差。

对,0001.0,01.0,1.0===εεε考虑一样问题。

解 (1)给出算法:为解b Ax =, 令U L D A --=, 其中][ij a A =, ),,,(2211nn a a a diag D = ,⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡------=-00001,21323121n n n n a a a a a a L,⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡------=-0000,122311312 n n n n a a a a a a U 利用Jacobi 迭代法, G-S 迭代法, SOR 迭代法解线性方程组, 均能够下步骤求解: step1给定初始向量x0=(0,0,...,0), 最大迭代次数N, 精度要求c, 令k=1 step2令x=B*x0+gstep3若||x-x0||2<c, 算法停止, 输出解和迭代次数k, 不然, 转step4step4若k>=N,算法停止, 迭代失败, 不然, 令x0=x, 转step2在Jacobi 迭代法中, B=D -1*(L+U),g=D -1*b在G-S 迭代法中, B=D -1*(L+U),g=D -1*b在SOR 迭代法中, B=(D-w*L)-1*[(1-w)*D+w*U],g=w*(D-w*L)-1*b另外, 在SOR 迭代法中, 上面算法step1中要给定松弛因子w, 其中0<w<2 为计算结果, 要求w=0.5。

《数值计算方法》上机实验报告

《数值计算方法》上机实验报告

《数值计算方法》上机实验报告华北电力大学实验名称数值il•算方法》上机实验课程名称数值计算方法专业班级:电力实08学生姓名:李超然学号:200801001008 成绩: 指导教师:郝育黔老师实验日期:2010年04月华北电力大学实验报告数值计算方法上机实验报吿一.各算法的算法原理及计算机程序框图1、牛顿法求解非线性方程*对于非线性方程,若已知根的一个近似值,将在处展开成一阶xxfx ()0, fx ()xkk泰勒公式"f 0 / 2 八八,fxfxfxxxxx 0 0 0 0 0 kkkk2!忽略高次项,有,fxfxfxxx 0 ()()(),,, kkk右端是直线方程,用这个直线方程来近似非线性方程。

将非线性方程的**根代入,即fx ()0, X ,* fxfxxx 0 0 0 0, ,, kkkfx 0 fx 0 0,解出fX 0 *k XX,, k' fx 0 k水将右端取为,则是比更接近于的近似值,即xxxxk, Ik, Ikfx ()k 八XX, Ikk* fx()k这就是牛顿迭代公式。

,2,计算机程序框图:,见,,3,输入变量、输出变量说明:X输入变量:迭代初值,迭代精度,迭代最大次数,\0输出变量:当前迭代次数,当前迭代值xkl,4,具体算例及求解结果:2/16华北电力大学实验报吿开始读入l>k/fx()0?,0fx 0 Oxx,,01* fx ()0XX,,,?10kk, ,1,kN, ?xx, 10输出迭代输出X输出奇异标志1失败标志,3,输入变量、输出变量说明: 结束例:导出计算的牛顿迭代公式,并il •算。

(课本P39例2-16) 115cc (0), 求解结果:10. 75000010.72383710. 72380510. 7238052、列主元素消去法求解线性方程组,1,算法原理:高斯消去法是利用现行方程组初等变换中的一种变换,即用一个不为零的数乘 -个 方程后加只另一个方程,使方程组变成同解的上三角方程组,然后再自下而上 对上三角3/16华北电力大学实验报告方程组求解。

大连理工复试内容及形式

大连理工复试内容及形式

大连理工大学学术型硕士研究生各学科、专业
复试内容及形式
院、系(部)名称:电子与信息工程学院 (7)
院、系(部)名称:电气工程及应用电子技术系 (10)
院、系(部)名称:应用数学系
院、系(部)名称:物理与光电工程学院
院、系(部)名称:运载工程与力学学部院、系(部)名称:工程力学系
院、系(部)名称:船舶工程学院
院、系(部)名称:汽车工程学院
院、系(部)名称:机械工程学院
院、系(部)名称:材料科学与工程学院
院、系(部)名称:土木水利学院
院、系(部)名称:化工学院
院、系(部)名称:电子与信息工程学院
院、系(部)名称:能源与动力学院
院、系(部)名称:人文社会科学学院
院、系(部)名称:电气工程及应用电子技术系
院、系(部)名称:外国语学院
院、系(部)名称:体育教育部
院、系(部)名称:建筑与艺术学院
院系名称:软件学院
院、系(部)名称:环境与生命学院
院、系(部)名称:马克思主义学院
院、系(部)名称:经济系。

线性代数 第二版 王希云 习题解答

线性代数 第二版 王希云 习题解答

解:由定义知,原式=0
2+ 2
× 3 + 0 × ( −1)
2+3
× ( −7 ) + 1× ( −1)
2+ 4
15 ×4 =
3. 计算 5 阶行列式
2 1 0 0 0 1 2 1 0 0 D5 = 0 1 2 1 0 0 0 1 2 1 0 0 0 1 2
原式
2 1 1 = 2 × ( −1) × 0 0
1 2 1 0
x+ y −y = −2 ( x 3 + y 3 ) y−x
a1 − b1 a1 − b2 a1 − bn a2 − b1 a2 − b2 a2 − bn = 0 an − b1 an − b2 an − bn
3. 由 n 阶行列式
1 1 1
1 1 1 1 =0 1 1
3. 写出四阶行列式 det( aij ) 中包含因子 a42 a23 的项,并指出正负号. 解;a11a23a42 a34 , − a4 a23a31a42 (全书解析请关注 v 信 宫中号 高校 课后习题) 4. 写出四阶行列式 det( aij ) 中所有取负号且包含因子 a23 的项. 解:取负号 N
( 0 + 0 + 0 +0 + n −1)
= ( −1)
n!
(5)原式 = 0 (不同行式不同列相乘, 每次均会出现 0, 故其和为 0) 6. 问
a11 0 0 a41
0 a22 a32 0
a14 a23 0 a11a22 a33a44 − a14 a23a32 a41 = a33 0 0 a44
行列式的值是不同行不同列之积的和对于后置来说对应和运算法则没有影响故同理可证性质五利用行列式的定义计算下列行列式

数值分析上机实验指导书

数值分析上机实验指导书

“数值计算方法”上机实验指导书实验一 误差分析实验1.1(病态问题)实验目的:算法有“优”与“劣”之分,问题也有“好”与“坏”之别。

对数值方法的研究而言,所谓坏问题就是问题本身对扰动敏感者,反之属于好问题。

通过本实验可获得一个初步体会。

数值分析的大部分研究课题中,如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。

病态问题要通过研究和构造特殊的算法来解决,当然一般要付出一些代价(如耗用更多的机器时间、占用更多的存储空间等)。

问题提出:考虑一个高次的代数多项式)1.1()()20()2)(1()(201∏=−=−−−=k k x x x x x p显然该多项式的全部根为1,2,…,20共计20个,且每个根都是单重的。

现考虑该多项式的一个扰动)2.1(0)(19=+x x p ε其中ε是一个非常小的数。

这相当于是对(1.1)中19x 的系数作一个小的扰动。

我们希望比较(1.1)和(1.2)根的差别,从而分析方程(1.1)的解对扰动的敏感性。

实验内容:为了实现方便,我们先介绍两个MATLAB 函数:“roots ”和“poly ”。

roots(a)u =其中若变量a 存储n+1维的向量,则该函数的输出u 为一个n 维的向量。

设a 的元素依次为121,,,+n a a a ,则输出u 的各分量是多项式方程01121=+++++−n n n n a x a x a x a的全部根;而函数 poly(v)b =的输出b 是一个n+1维向量,它是以n 维向量v 的各分量为根的多项式的系数(从高到低排列)。

可见“roots ”和“poly ”是两个互逆的运算函数。

))20:1((;)2();21,1(;000000001.0ve poly roots ess ve zeros ve ess +===上述简单的MATLAB 程序便得到(1.2)的全部根,程序中的“ess ”即是(1.2)中的ε。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

(2)估计5到20阶Hilbert 矩阵的∞范数条件数(3)设n n R A ⨯∈⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡------=111111111011001,先随机地选取n R x ∈,并计算出x A b n =;然后再用列主元Gauss 消去法求解该方程组,假定计算解为∧x 。

试对n 从5到30估计计算解∧x 的精度,并且与真实相对误差作比较。

解(1)分析:利用for 使n 从5循环到20,利用()hilb 函数得到Hilbert 矩阵A ;先将算法编制成通用的子程序,利用算法编成的子程序)(B opt v =,对TA B -=求解,得到∞-1A的一个估计值v v =~;再利用inf),(A norm 得到∞A ;则条件数inf),(1A norm v A A K *==∞∞-。

另,矩阵A 的∞范数条件数可由inf),(A cond 直接算出,两者可进行比较。

程序为1 算法编成的子程序)(B opt v =function v=opt(B)k=1;;n=length(B);x=1./n*ones(n,1); while k==1 w=B*x;v=sign(w); z=B'*v;if norm(z,inf)<=z'*x v=norm(w,1);…k=0; elsex=zeros(n,1);[s,t]=max(abs(z)); x(t)=1; k=1; end end2 问题(1)求解 ex2_1for n=5:20A=hilb(n);B=inv(A.');v=opt(B);K1=v*norm(A,inf);(K2=cond(A,inf);disp(['n=',num2str(n)])disp(['估计条件数为',num2str(K1)])disp(['实际条件数为',num2str(K2)])end计算结果为n=5估计条件数为943656—实际条件数为943656n=6估计条件数为.0028实际条件数为.0028n=7估计条件数为实际条件数为n=8|估计条件数为实际条件数为n=9估计条件数为.422实际条件数为.422n=10估计条件数为实际条件数为,n=11估计条件数为344实际条件数为344Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .In ex2_1 at 6n=12估计条件数为+16实际条件数为+16Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .> In ex2_1 at 3^Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .> In cond at 47In ex2_1 at 6n=13估计条件数为+18实际条件数为+18Warning: Matrix is close to singular or badly scaled. Results may be ;inaccurate. RCOND = .> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .> In cond at 47In ex2_1 at 6n=14估计条件数为+17、实际条件数为+17Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .> In cond at 47In ex2_1 at 6。

n=15估计条件数为+17实际条件数为+17Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .In ex2_1 at 6n=16估计条件数为+17实际条件数为+17Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .> In ex2_1 at 3&Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .> In cond at 47In ex2_1 at 6n=17估计条件数为+18实际条件数为+18Warning: Matrix is close to singular or badly scaled. Results may be !inaccurate. RCOND = .> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .> In cond at 47In ex2_1 at 6n=18估计条件数为+18【实际条件数为+18Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .> In cond at 47In ex2_1 at 6}n=19估计条件数为+18实际条件数为+18Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .In ex2_1 at 6 n=20估计条件数为+18 实际条件数为+18结果分析随着矩阵阶数增加,估计值误差开始出现,20,17,16,15=n 时估计条件数与实际值存在误差;且条件数很大,Hilbert 矩阵为病态的。

解(2)分析:先根据题目要求,利用for 和()rand 使n 从5循环到30,作出A 和随机的x ,并计算出Ax b =;然后再利用第一章习题中得到的)(],,[A GaussCol P U L =和),,,,(P U L b A Gauss x =用列主元Gauss 消去法求解该方程组,假定计算解为1x ,得1*x A b r -=,利用第(1)问所得函数)).'((A inv opt 计算∞-1A 的一个估计值,利用inf),(*norm 计算A b r ,,的无穷范数,则1x 的相对误差估计为))/norm(b,(A,A.'))*norm )*opt(inv(norm(r,p inf inf inf 1=,真实相对误差为))/norm(x,,norm(x-x p inf inf 12=。

\程序为1 列主元Gauss 消去法求解该方程组的程序为 A 的LU 分解:function [L,U,P]=GaussCol(A) n=length(A); for k=1:n-1[s,t]=max(abs(A(k:n,k))); p=t+k-1; $temp=A(k,1:n);A(k,1:n)=A(p,1:n); A(p,1:n)=temp;u(k)=p;if A(k,k)~=0A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);@elsebreak;endendL=tril(A);U=triu(A);L=L-diag(diag(L))+diag(ones(1,n));P=eye(n);…for i=1:n-1temp=P(i,:);P(i,:)=P(u(i),:);P(u(i),:)=temp;endend高斯消去法解线性方程组?function x=Gauss(A,b,L,U,P)if nargin<5P=eye(length(A));endn=length(A);b=P*b;for j=1:n-1b(j)=b(j)/L(j,j);}b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endb(n)=b(n)/L(n,n);y=b;for j=n:-1:2y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);'endy(1)=y(1)/U(1,1);x=y;end2 问题(2)求解ex2_2for n=5:30A=2*eye(n)+tril(-1*ones(n)); A(1:n-1,n)=ones(n-1,1); x=100*rand(n,1);b=A*x;[L,U,P]=GaussCol(A); x1=Gauss(A,b,L,U,P);r=b-A*x1;p1=norm(r,inf)*opt(inv(A.'))*norm(A,inf)/norm(b,inf); p2=norm(x-x1,inf)/norm(x,inf);disp(['n=',num2str(n)])disp(['估计相对误差为',num2str(p1)])disp(['实际相对误差为',num2str(p2)]);y1(n-4)=p1;y2(n-4)=p2;endplot(5:30,y1,5:30,y2)legend('估计相对误差','实际相对误差')计算结果为n=5估计相对误差为实际相对误差为.n=6估计相对误差为实际相对误差为n=7估计相对误差为实际相对误差为n=8估计相对误差为^实际相对误差为n=9估计相对误差为实际相对误差为n=10估计相对误差为实际相对误差为n=11(估计相对误差为实际相对误差为n=12估计相对误差为实际相对误差为n=13估计相对误差为实际相对误差为`n=14估计相对误差为实际相对误差为n=15估计相对误差为实际相对误差为n=16估计相对误差为—实际相对误差为n=17估计相对误差为实际相对误差为n=18估计相对误差为实际相对误差为n=19;估计相对误差为实际相对误差为n=20估计相对误差为实际相对误差为n=21估计相对误差为实际相对误差为n=22估计相对误差为实际相对误差为n=23估计相对误差为实际相对误差为n=24估计相对误差为实际相对误差为n=25估计相对误差为实际相对误差为n=26估计相对误差为实际相对误差为n=27-8估计相对误差为实际相对误差为n=28估计相对误差为实际相对误差为n=29估计相对误差为实际相对误差为n=30估计相对误差为实际相对误差为结果分析n较小时估计的较好,随着n的增大估计值误差增大。

相关文档
最新文档