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

合集下载

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

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

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

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

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

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

程序为1 算法2、5、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 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估计条件数为29070279、0028实际条件数为29070279、0028n=7估计条件数为985194887、5079实际条件数为985194887、5079n=8估计条件数为33872789099、7717实际条件数为33872789099、7717n=9估计条件数为16、422实际条件数为16、422n=10估计条件数为35353368771750、67实际条件数为35353368771750、67n=11估计条件数为1232433965549344实际条件数为1232433965549344Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、547634e-17、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、547634e-17、> In cond at 47In ex2_1 at 6n=12估计条件数为3、9245e+16实际条件数为3、9245e+16Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 7、847381e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 7、847381e-19、> In cond at 47In ex2_1 at 6n=13估计条件数为1、2727e+18实际条件数为1、2727e+18Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、246123e-18、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、246123e-18、> In cond at 47In ex2_1 at 6n=14估计条件数为4、8374e+17实际条件数为4、8374e+17Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 8、491876e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 8、491876e-19、> In cond at 47In ex2_1 at 6n=15估计条件数为4、6331e+17实际条件数为5、234289848563619e+17Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 9、137489e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 9、137489e-19、> In cond at 47In ex2_1 at 6n=16估计条件数为8、3166e+17实际条件数为8、3167e+17Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 6、244518e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 6、244518e-19、 > In cond at 47 In ex2_1 at 6 n=17估计条件数为1、43e+18 实际条件数为1、43e+18Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、693737e-19、 > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、693737e-19、 > In cond at 47 In ex2_1 at 6 n=18估计条件数为2、5551e+18 实际条件数为2、8893e+18Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、264685e-19、 > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、264685e-19、 > In cond at 47 In ex2_1 at 6 n=19估计条件数为2、411858563109357e+18 实际条件数为2、411858563109357e+18Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 1、351364e-19、 > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 1、351364e-19、 > In cond at 47 In ex2_1 at 6 n=20估计条件数为2、31633670586674e+18 实际条件数为6、37335273308473e+18结果分析随着矩阵阶数增加,估计值误差开始出现,20,17,16,15 n 时估计条件数与实际值存在误差;且条件数很大,Hilbert 矩阵为病态的。

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

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

数值线性代数第二版徐树方高立张平文上机习题第一章实验报告(供参考)上机习题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代码%算法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%算法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); else break;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;end。

数值代数上机实验报告

数值代数上机实验报告

数值代数上机实验报告试验项目名称:平方根法与改进平方根法实验内容:先用你熟悉的计算机语言将平方根法和改进平方根法编写成通用的子程序,然后用你编写的程序求解对称正定方程组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<>。

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华北电力大学实验报告方程组求解。

数值线性代数(徐树芳老师)第一章答案

数值线性代数(徐树芳老师)第一章答案

数值线性代数(徐树芳⽼师)第⼀章答案习题11.求下三⾓阵的逆矩阵的详细算法。

[解] 设下三⾓矩阵L的逆矩阵为T我们可以使⽤待定法,求出矩阵T的各列向量。

为此我们将T按列分块如下:注意到我们只需运⽤算法1·1·1,逐⼀求解⽅程便可求得[注意]考虑到内存空间的节省,我们可以置结果矩阵T的初始状态为单位矩阵。

这样,我们便得到如下具体的算法:算法(求解下三⾓矩阵L的逆矩阵T,前代法)2.设为两个上三⾓矩阵,⽽且线性⽅程组是⾮奇异的,试给出⼀种运算量为的算法,求解该⽅程组。

[解]因,故为求解线性⽅程组,可先求得上三⾓矩阵T的逆矩阵,依照上题的思想我们很容易得到计算的算法。

于是对该问题我们有如下解题的步骤:(1)计算上三⾓矩阵T的逆矩阵,算法如下:算法 1(求解上三⾓矩阵的逆矩阵,回代法。

该算法的的运算量为)(2)计算上三⾓矩阵。

运算量⼤约为.(3)⽤回代法求解⽅程组:.运算量为;(4)⽤回代法求解⽅程组:运算量为。

算法总运算量⼤约为:3.证明:如果是⼀个Gauss变换,则也是⼀个Gauss 变换。

[解]按Gauss变换矩阵的定义,易知矩阵是Gauss变换。

下⾯我们只需证明它是Gauss变换的逆矩阵。

事实上注意到,则显然有从⽽有4.确定⼀个Gauss变换L,使量的第⼆⾏加上第⼀⾏的2倍;使向量的第三⾏加上第⼀⾏的2倍。

于是Gauss变换如下5.证明:如果有三⾓分解,并且是⾮奇异的,那么定理1·1·2中的L和U都是唯⼀的。

[证明]设,其中都是单位下三⾓阵,都是上三⾓阵。

因为A⾮奇异的,于是注意到,单位下三⾓阵的逆仍是单位下三⾓阵,两个单位下三⾓阵的乘积仍是单位下三⾓阵;上三⾓阵的逆仍是上三⾓阵,两个上三⾓阵的乘积仍是上三⾓阵。

因此,上述等将是⼀个单位下三⾓阵与⼀个上三⾓阵相等,故此,它们都必是单位矩阵。

即,从⽽即A的LU分解是唯⼀的。

6.设的定义如下证明A有满⾜的三⾓分解。

线性代数第二版第一章1

线性代数第二版第一章1

a11 a12 = a11a22 - a12a21 a21 a22
并称之为二阶行列式.
备注
定义
由四个数排成二行二列(横排称行、 由四个数排成二行二列(横排称行、竖排
称列) 称列)的数表
a11 a12 a21 a22 ( 4)
表达式 a11a22 − a12 a21称为数表( 4)所确定的二阶 称为数表( a11 行列式, 行列式,并记作 a21
n
= l
1
l
2
L l
n
,
例5 计算对角行列式
0 0 0 4 0 0 3 0 0 2 0 0 1 0 0 0
五.小结
(1) 上三角形行列式 (主对角线下侧元素都为0) 主对角线下侧元素都为 )
a11 a12 L a1n a2 n = a11a22 L ann M ann
0 a22 L D= M M O 0 0 L
其中每一条实线上的三个元素的乘积带正号, 其中每一条实线上的三个元素的乘积带正号,每一 条虚线上的三个元素的乘积带负号, 条虚线上的三个元素的乘积带负号,所得六项的代 数和就是三阶行列式的展开式. 数和就是三阶行列式的展开式.
注意
例 2
计算三阶行列式
1 D= 2 - 1
- 2 1 1
1 - 3 - 1
a11 a21 a31
a12 a22 a32
a13 a23 = a11a22 a33 + a12 a23 a31 + a13 a21a32 a33 - a11a23 a32 - a12 a21a33 - a13 a22 a31
= a11 (a22 a33 - a23 a32 ) - a12 (a21a33 - a23 a31 ) + a13 (a21a32 - a22 a31 )

代数学引论高教第二版答案(第一章)

代数学引论高教第二版答案(第一章)

1. 如果群 G 中,对任意元素 a,b 有(ab)2=a2b2,则 G 为交换群. 证明: 对任意 a,b G,由结合律我们可得到(ab)2=a(ba)b, a2b2=a(ab)b 再由已知条件以及消去律得到ba=ab, 由此可见群 G 为交换群.2. 如果群 G 中,每个元素 a 都适合 a2=e, 则 G 为交换群. 证明: [方法 1] 对任意 a,b G,ba=bae=ba(ab)2=ba(ab)(ab) =ba2b(ab)=beb(ab)=b2(ab)=e(ab)=ab 因此 G 为交换群. [方法 2] 对任意 a,b G,a2b2=e=(ab)2, 由上一题的结论可知 G 为交换群.3. 设 G 是一非空的有限集合,其中定义了一个乘法 ab,适合条件: (1) a(bc)=(ab)c; (2) 由 ab=ac 推出 b=c; (3) 由 ac=bc 推出 a=b;证明 G 在该乘法下成一群. 证明:[方法 1]设 G={a1,a2,…,an},k 是 1,2,…,n 中某一个数字,由(2)可知若 i j(I,j=1,2,…,n),有 akai ak aj------------<1> aiak aj ak------------<2>再由乘法的封闭性可知 G={a1,a2,…,an}={aka1, aka2,…, akan}------------<3> G={a1,a2,…,an}={a1ak, a2ak,…, anak}------------<4>由<1>和<3>知对任意 at G, 存在 am G,使得 akam=at.由<2>和<4>知对任意 at G, 存在 as G,使得 asak=at.由下一题的结论可知 G 在该乘法下成一群.1下面用另一种方法证明,这种方法看起来有些长但思路比较清楚。

线性代数简明教程-第二版-答案

线性代数简明教程-第二版-答案

3. 求下列排列的逆序数
(1) (315624) 6
(2) (13(2n 1)24(2n)) n(n 1)
2
4. 计算下列行列式
2500 350
55
(1)
500 70
1500 70
31
35000(5 15) 350000
a11 a12 a21 a22 (2) 0 0
00 00 a33 a34
2 0
0
0 1 0
0
0 1
4
4 0 0
0 3 0
0 0 2
2 0 0
0 3 0
0
0 1
2
12.设
1 2 3
A
0 3 0
2 2 1
2 0 2
利用初等行变换求 A1
2
1 1 1
,
1 2 3 2 1 0 0 0
(
A
E)
0 3 0
2 2 1
2 0 2
1 0 1 0 0
5、已知两个线性变换
x2
x1 2 y1 y3 2 y1 3y2
2 y3
x3 4 y1 y2 5 y3
y1 3z1 z2 y2 2z1 z3
,
y3 z2 3z3
求从 z1, z2 , z3 到 x1, x2 , x3 的 线性变换
分析:X AY ,Y BZ ,
1 1
0 0
0 0
1 0
10
1 2 3 2 1 0 0 0
~r3 3r1 0
0 0
2 4 1
2 9 2
1 0 1 0 0
5 1
3 0
0 0
1 0
0 1
1 2 3 2 1 0 0 0

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

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

第四章上机习题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 ,否则,转step4 step4若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。

数值分析上机实验指导书

数值分析上机实验指导书

“数值计算方法”上机实验指导书实验一 误差分析实验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)中的ε。

线性代数简明教程,(第二版)科学出版社课后答案第一章矩阵习题答案

线性代数简明教程,(第二版)科学出版社课后答案第一章矩阵习题答案

0

0 0
1 0 0
0 0 0 1 0 1
1 0
0 1
3 6
1 1
1 2
610
1 0 0 0 2 1 1 4


~r1 2r2
0

0 0
1 0 0
0 1 0
0 0 1
0 3 6
1 1 1
0 1 2
1
6 10

(E A1 )
1 2 3 1
r2 2r1
r3 r1
1 2 3 1 0 0 0 0
0 0 0 0
3 1 4 2 2
(2) B


1 1 1
0 2 4
1 1 13 3 3
0
4 0

1 0 1 1 0

X X
11 21
X12 X 22


E


En O
O Es

AX21 En

A X 22 BX11

O O
BX12 Es

X 21 A1 X 22 O
X11 O
X12 B1
方法二:

Ons Bs
AB A 2B
B1 ( A 2E)1 A
方法一:先求 ( A 2E)1, 再求 ( A 2E)1 A
方法二:直接求 ( A 2E)1 A

(( A 2E)A)
(E( A 2E)1 A)
其中
1 0 0
(A

2E)1


2 0

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

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

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

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

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

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

程序为1 算法2、5、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 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估计条件数为29070279、0028实际条件数为29070279、0028n=7估计条件数为985194887、5079实际条件数为985194887、5079n=8估计条件数为33872789099、7717实际条件数为33872789099、7717n=9估计条件数为16、422实际条件数为16、422n=10估计条件数为35353368771750、67实际条件数为35353368771750、67n=11估计条件数为1232433965549344实际条件数为1232433965549344Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、547634e-17、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、547634e-17、> In cond at 47In ex2_1 at 6n=12估计条件数为3、9245e+16实际条件数为3、9245e+16Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 7、847381e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 7、847381e-19、> In cond at 47In ex2_1 at 6n=13估计条件数为1、2727e+18实际条件数为1、2727e+18Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、246123e-18、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、246123e-18、> In cond at 47In ex2_1 at 6n=14估计条件数为4、8374e+17实际条件数为4、8374e+17Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 8、491876e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 8、491876e-19、> In cond at 47In ex2_1 at 6n=15估计条件数为4、6331e+17实际条件数为5、234289848563619e+17Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 9、137489e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 9、137489e-19、> In cond at 47In ex2_1 at 6n=16估计条件数为8、3166e+17实际条件数为8、3167e+17Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 6、244518e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 6、244518e-19、 > In cond at 47 In ex2_1 at 6 n=17估计条件数为1、43e+18 实际条件数为1、43e+18Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、693737e-19、 > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、693737e-19、 > In cond at 47 In ex2_1 at 6 n=18估计条件数为2、5551e+18 实际条件数为2、8893e+18Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、264685e-19、 > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、264685e-19、 > In cond at 47 In ex2_1 at 6 n=19估计条件数为2、411858563109357e+18 实际条件数为2、411858563109357e+18Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 1、351364e-19、 > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 1、351364e-19、 > In cond at 47 In ex2_1 at 6 n=20估计条件数为2、31633670586674e+18 实际条件数为6、37335273308473e+18结果分析随着矩阵阶数增加,估计值误差开始出现,20,17,16,15 n 时估计条件数与实际值存在误差;且条件数很大,Hilbert 矩阵为病态的。

线性代数 高教第二版 卢刚 第一章 课后答案

线性代数 高教第二版 卢刚 第一章 课后答案

7

第一章《线性代数》习题解答
10 2 3 4
10 2 3 4
(6)原
第2,3,4行都加到第1行
10
3
4
1 0 第2,3,4行都减去第1行
10 4 1 2
0
1 1 -3 2 -2 -2
10 1 2 3
0 -1 -1 -1
1 1 3 10 2 2 2 160.
按第1列展开
14. A2 A [1 (B I )]2 1 (B I ) 1 (B 2 2BI I 2 ) 1 (B I )
2
2
4
2
(B 2 2B I ) 2(B I ) B 2 I.
n
n
n
15.(1) tr( A B) (aii bii ) aii bii tr( A) tr(B) .
第2列减去第1列 28092 1000
28092 1
5 -1 3 5 -1 3
第2行的( 98) 倍加到第3行
按第1行展开
(4) 原
2 2 2=2 1 1 1 = 8 .
0 73 0 73
(5)利用 P22 例 6 的结论.原 [5 (4 1) 1](5 1)3 512 .
第 7 页 共 20 页
0 1 0 0
0 0 1 0
0 0 0 1
(4)(直接计算即可)令
A
0 0 0
0 0 0
1 0 0
0 10
,则
A2
0
0 0
0 0 0
0 0 0
1
0 0

A3
0
0 0
0 0 0
0 0 0
0
0 0

0 0 0 0

数值线性代数大作业报告

数值线性代数大作业报告

数值线性代数实验大报告指导老师:赵国忠姓名:1108300001 刘帅1108300004 王敏1108300032 郭蒙一、实验名称:16题P75上机习题二、实验目的:编制通用的子程序,完成习题的计算任务三、实验内容与要求:P75上机习题先用熟悉的计算机语言将算法2.5.1编制成通用的子程序,然后再用所编制的子程序完成下面两个计算任务:(1) 估计5到20阶Hilbert 矩阵的无穷范数条件数。

(2) 设A n = 11...111................1-1 (01)-- 先随机地选取x ∈R n ,并计算出b=An x;然后再用列主元Gauss 消去法求解该方程组,假定计算解为∧x .试对n 从5到30估计计算解∧x 的精度,并且与真实的相对误差作比较。

四、 实验原理:(1)矩阵范数(martix norm )是数学上向量范数对矩阵的一个自然推广。

利用for循环和cond (a )Hilbert 求解Hilbert 矩阵的无穷范数,再利用norm(a,inf)求矩阵的无穷范数条件数。

(2)本题分为4步来求解。

先运用rand 随机选取x ∈R n,输入A n 矩阵,编制一个M 文件计算出b 。

第二步用列主元高斯消去法求解出方程的解X2。

第三步建立M 文件: soluerr.m 估计计算解∧x 的精度。

第四步, 建立M 文件: bijiao.m ,与真实相对误差作比较。

五、 实验过程:(1)程序:clearfor n=5:20for i=1:nfor j=1:na(i,j)=1/(i+j-1);endendc=cond(a);f=norm(c,inf);fprintf('n=%3.0f\nnorm(c,inf)%e\n',n,f) end运行结果:n= 5norm(c,inf)4.766073e+005n= 6norm(c,inf)1.495106e+007n= 7norm(c,inf)4.753674e+008n= 8norm(c,inf)1.525758e+010n= 9norm(c,inf)4.931542e+011n= 10norm(c,inf)1.602467e+013n= 11norm(c,inf)5.224376e+014n= 12norm(c,inf)1.698855e+016n= 13norm(c,inf)3.459404e+017n= 14norm(c,inf)4.696757e+017n= 15norm(c,inf)2.569881e+017n= 16norm(c,inf)7.356249e+017n= 17norm(c,inf)4.362844e+017n= 18norm(c,inf)1.229633e+018n= 19norm(c,inf)9.759023e+017n= 20norm(c,inf)1.644051e+018(2)程序:M文件:matrix1.mfunction [a,b,x1]=matrix1(n) format longA1=-1*ones(n,n)A2=tril(A1)for i=1:nA2(i,i)=1endA2(:,n)=1a=A2x1=rand(n,1)b=A2*x1end运行结果:>> A1 =-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1A2 =-1 0 0 0 0 -1 -1 0 0 0 -1 -1 -1 0 0 -1 -1 -1 -1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0 -1 -1 0 0 0 -1 -1 -1 0 0 -1 -1 -1 -1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0-1 -1 -1 0 0 -1 -1 -1 -1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0 -1 1 0 0 0 -1 -1 1 0 0 -1 -1 -1 -1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0 -1 1 0 0 0 -1 -1 1 0 0 -1 -1 -1 1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0 -1 1 0 0 0 -1 -1 1 0 0 -1 -1 -1 1 0 -1 -1 -1 -1 1A2 =1 0 0 0 1 -1 1 0 0 1 -1 -1 1 0 1 -1 -1 -1 1 1 -1 -1 -1 -1 1a =-1 1 0 0 1 -1 -1 1 0 1 -1 -1 -1 1 1 -1 -1 -1 -1 1x1 =0.8147236863931790.9057919370756190.1269868162935060.9133758561390190.632359246225410b =1.4470829326185890.723427496907850-0.961169560949882-0.301767337397875-2.128519049675914a =1 0 0 0 1 -1 1 0 0 1 -1 -1 1 0 1 -1 -1 -1 1 1 -1 -1 -1 -1 1b =1.4470829326185890.723427496907850-0.961169560949882-0.301767337397875-2.128519049675914x1 =0.8147236863931790.9057919370756190.1269868162935060.9133758561390190.632359246225410M文件:LZYgauss.mfunction[x2]=LZYgauss(a,b)format longn=length(a);x2=zeros(n,1);a=[a b];for k=1:n-1max=k;for i=k+1:nif a(i,k)>a(max,k)max=i;endendtemp=a(k,k:n+1);a(k,k:n+1)=a(max,k:n+1);a(max,k:n+1)=temp;for i=k+1:na(i,k)=-a(i,k)/a(k,k);a(i,k+1:n+1)=a(i,k+1:n+1)+a(i,k)*a(k,k+1:n+1);endendx2(n,1)=a(n,n+1)/a(n,n);for i=n-1:-1:1sum=0;for j=i+1:nsum=sum+x2(j,1)*a(i,j);endx2(i,1)=(a(i,n+1)-sum)/a(i,i);end运行结果:>> LZYgauss(a,b)ans =0.8147236863931790.9057919370756190.1269868162935060.913375856139020 0.632359246225410估计计算解x的精度:M文件: soluerr.mfunction [x,error]=soluerr(a,b)format long%估计计算解的精度% 算法:列主元Gauss消去法,其中% A --- 系数矩阵% b-右端项% index --- index=0表示计算成败;index=1表示计算成功%输出结果:error--本算法给出的计算解的估计% normA--逆矩阵无穷范数估计% rnorm--计算解的残量[n,m]=size(a); nb=length(b);if n~=merror('The rows and columns of matrix a must be equal!');return;endif m~=nberror('The columns of a must be equal the dimension of b!'); return;endindex=1;%列主元矩阵三角分解[L,U,u,index_col]=Gauss_col(a);%解下三角方程组Ly=Pb[y,index_low]=Gauss_low(L,b(u));%解上三角方程组Ux=y[x,index_upp]=Gauss_upp(U,y);%输出数值解xpause(0.3)%估计矩阵逆的无穷大范数normA=normAinv(L,U,u);%估计计算解的残量rnorm=norm(b-a*x,inf);%计算右端项bnorm=norm(b,inf);%计算矩阵A的范数Anorm=norm(a,inf);%计算解的精度error=normA*Anorm*rnorm/bnorm;运行结果:>> [x,error]=soluerr(a,b)x =0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410x =0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410error =5.215941218821592e-016X1与真实相对误差做比较M文件:bijiao.mfunction [x,x1,error,error1]=bijiao(a,b,n) [a,x1]=matrix1(n)[x,error]=soluerr(a,b)error1=abs((x-x1)/x1)运行结果:>>x =0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410x1 =0.9578935505663481.6132601689718680.629241553693582-0.799716694069468-1.770467991515150error =5.215941218821592e-016error1 =Columns 1 through 40 0 0 00 0 0 00 0 0 00 0 0 00 0 0 0Column 50.0808655478999350.3995939126190040.2836847318376260.9675930648949141.357170674226221六、结果分析:1、矩阵范数(martix norm)是数学上向量范数对矩阵的一个自然推广。

数值线性代数实验

数值线性代数实验

数值线性代数实验题目:数值线性代数专业:信息与计算科学班级:班姓名:山东科技大学2013年 1 月16日实验报告说明学院:信息学院专业:信息班级10-2 姓名:一、主要参考资料:(1)《Matlab数值计算-案例分析》北京航空出版(2)《Matlab数值分析》机械工业出版二、课程设计应解决的主要问题:(1)平方根(2)QR方法(3)最小二乘法三、应用软件:(1)Matlab7.0(2)数学公式编辑器四、发出日期:课程设计完成日期:指导教师签字:系主任签字:指导教师对课程设计的评语指导教师签字:年月日一、问题描述先用你所熟悉的计算机语言将平方根和改进的平方根法编成写通用的子程序,然后用你编写的程序求解对称正定方程组b x =A ,其中 (1)b 随机的选取,系数矩阵位100阶矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡1011101110111011101110(2)系数矩阵为40阶Hilbert 矩阵,即系数矩阵A 的第i 行第j 列元素为11-+=j i a ij ,向量b 的第i 个分量为∑=-+=nj i j i b 111。

二、分析与程序1. 平方根法函数程序如下: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 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=bL=zeros(n,n);D=diag(n,0);S=L*D;for i=1:n %L的主对角元素均为1L(i,i)=1;endfor i=1:nfor j=1:nif (eig(A)<=0)|(A(i,j)~=A(j,i))disp('wrong');break;endendendD(1,1)=A(1,1);for 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=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); endfor i=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n));end2.改进平方根法函数程序如下: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);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=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'); end3.调用函数解题: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[x,b]=pingfanggenfa(A,b)b=gaijinpinfanggenfa(A,b)4.运行结果:请输入矩阵维数: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.000012.0000 -69.000068.000021.000017.0000 -50.7188-8.7500-8.0000 112.00006.0000 -68.750022.000044.0000 -28.00008.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.3285-0.44380.4621-0.25130.05650.0000-0.00510.0071-0.0027-0.00310.0036-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.00120.0139-0.09540.4208-1.21012.0624-1.0394-3.33436.2567-0.2463-7.45942.80303.69900.7277-1.7484-0.4854-3.60100.25325.1862-2.12991.44100.8738-4.56541.04224.0920-2.7764-2.2148-0.89530.36654.89671.04160.1281 -4.3387 -1.1902 -2.8334 8.4610 -3.6008一、问题描述先用你所熟悉的计算机语言将算法2.5.1编成写通用的子程序,然后用你编写的程序完成下面两个计算任务:(1) 估计5到20阶Hilbert 矩阵的∞范数条件数;(2)设n *n n R 11-1-1-111-1-101-1001A ∈⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=先随机的选取n R x ∈,并计算出x A b n =;然后再用列主元Gauss 消去法求解该方程组,假定计算解为xˆ。

[VIP专享]数值线性代数(徐树芳老师)第一章答案

[VIP专享]数值线性代数(徐树芳老师)第一章答案

习题11.求下三角阵的逆矩阵的详细算法。

[解] 设下三角矩阵L的逆矩阵为T我们可以使用待定法,求出矩阵T的各列向量。

为此我们将T按列分块如下:注意到我们只需运用算法1·1·1,逐一求解方程便可求得[注意]考虑到内存空间的节省,我们可以置结果矩阵T的初始状态为单位矩阵。

这样,我们便得到如下具体的算法:算法(求解下三角矩阵L的逆矩阵T,前代法)2.设为两个上三角矩阵,而且线性方程组是非奇异的,试给出一种运算量为的算法,求解该方程组。

[解]因,故为求解线性方程组,可先求得上三角矩阵T的逆矩阵,依照上题的思想我们很容易得到计算的算法。

于是对该问题我们有如下解题的步骤:(1)计算上三角矩阵T的逆矩阵,算法如下:算法 1(求解上三角矩阵的逆矩阵,回代法。

该算法的的运算量为)(2)计算上三角矩阵。

运算量大约为.(3)用回代法求解方程组:.运算量为;(4)用回代法求解方程组:运算量为。

算法总运算量大约为:3.证明:如果是一个Gauss变换,则也是一个Gauss变换。

[解]按Gauss变换矩阵的定义,易知矩阵是Gauss变换。

下面我们只需证明它是Gauss变换的逆矩阵。

事实上注意到,则显然有从而有4.确定一个Gauss变换L,使[解] 比较比较向量和可以发现Gauss变换L应具有功能:使向量的第二行加上第一行的2倍;使向量的第三行加上第一行的2倍。

于是Gauss变换如下5.证明:如果有三角分解,并且是非奇异的,那么定理1·1·2中的L和U都是唯一的。

[证明]设,其中都是单位下三角阵,都是上三角阵。

因为A非奇异的,于是注意到,单位下三角阵的逆仍是单位下三角阵,两个单位下三角阵的乘积仍是单位下三角阵;上三角阵的逆仍是上三角阵,两个上三角阵的乘积仍是上三角阵。

因此,上述等将是一个单位下三角阵与一个上三角阵相等,故此,它们都必是单位矩阵。

即,从而即A的LU分解是唯一的。

6.设的定义如下证明A有满足的三角分解。

数值代数上机实验报告

数值代数上机实验报告

数值代数课程设计实验报告姓名: 班级: 学号: 实验日期:一、实验名称 代数的数值解法 二、实验环境实验一、平方根法与改良平方根法一、实验要求:用熟悉的运算机语言将不选主元和列主元Gasuss 消元法编写成通用的子程序,然后用编写的程序求解以下方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--⨯14151515157681681681681681612321n n n n n x x x x x x 用所编的程序别离求解40、84、120阶方程组的解。

二、算法描述及实验步骤GAuss 程序如下:(1)求A 的三角分解:LU A =;(2)求解b y =L 得y ; (3)求解y x =U 得x ;列主元Gasuss 消元法程序如下: 1求A 的列主元分解:LU PA =;2求解b y P L =得y ;3求解y x U 得x ;三、调试进程及实验结果:%----------------方程系数---------------->> A1=Sanduijiaozhen(8,6,1,40); >> A2=Sanduijiaozhen(8,6,1,84); >> A3=Sanduijiaozhen(8,6,1,120); >> b1(1)=7;b2(1)=7;b3(1)=7; >> for i=2:39 b1(i)=15; end>> b1(40)=14; >> for i=2:83 b2(i)=15; end>> b2(40)=14; >> for i=2:119 b1(i)=15; end>> b3(120)=14;%----------------方程解---------------->> x11=GAuss(A1,b1') >> x12=GAuss Zhu (A1,b1') >> x21=GAuss(A2,b2') >> x22=GAuss Zhu (A3,b3') >> x31=GAuss(A3,b3') >> x32=GAuss Zhu (A3,b3') 运行结果:(n=40) GAuss 消元法的解即为x11 =列主元GAuss消元法的解即为x12 =1 1 1 1 1 1 1 1 1 1 111111111111111111111111111111六、源程序:function A=Sanduijiaozhen(a,b,c,n)%生成n阶以a,b,c为元素的三对角阵A=diag(b*ones(1,n),0)+diag(c*ones(1,n-1),1)+diag(a*ones(1,n-1),-1);function x=GAuss(A,b)n=length(b);x=b;%-------分解---------------for i=1:n-1for j=i+1:nmi=A(j,i)/A(i,i);b(j)=b(j)-mi*b(i);for k=i:nA(j,k)=A(j,k)-mi*A(i,k);endAB=[A,b];endend%-----------回代------------------x(n)=b(n)/A(n,n);for i=n-1:-1:1s=0;for j=i+1:ns=s+A(i,j)*x(j);endx(i)=(b(i)-s)/A(i,i);endfunction x=GAussZhu(A,b)n=length(b);x=b;%----------------------选主元---------------------for k=1:n-1a_max=0;for i=k:nif abs(A(i,k))>a_maxa_max=abs(A(i,k));r=i;endendif r>kfor j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%--------------消元-----------------for i=k+1:nm=A(i,k)/A(k,k);for j=k:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);endendif abs(A(n,n))==0return;endAbZhu=[A,b];%----------------回代-----------------------x(n)=b(n)/A(n,n);for i=n-1:-1:1for j=i+1:nb(i)=b(i)-A(i,j)*x(j);endx(i)=b(i)/A(i,i);end实验二、平方根法与改良平方根法一、实验要求:用运算机语言将平方根法和改良的平方根法编成通用的子程序,然后用编写的程序求解对称正定方程组100阶方程组AX=b,二、算法描述及实验步骤:平方根法函数程序如下:一、求A 的Cholesky 分解:L L A T=; 二、求解b y =L 得y ;3、求解y x =TL 得x ;改良平方根法函数程序如下:一、求A 的Cholesky 分解:T=LDL A ; 二、求解b y =L 得y ; 3、求解y x =TDL 得x ;三、调试进程及实验结果:clear;clc;%----------------方程系数---------------->> A=Sanduijiaozhen(1,10,1,100); >> b(1)=11; >> for i=2:99 b(i)=12; end>> b(100)=11;>> x1=Cholesky(A,b') >> x2=GJCholesky(A,b')运行结果:平方根法的解即为 x1 =改良平方根法解得的解即为x2 =四、源程序:function x=Cholesky(A,b)n=size(A);n=n(1);% x=A^-1*b;% disp('Matlab自带解即为x');%-----------------Cholesky分解-------------------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%------------------前代法求解Ly=b----------------------------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);%-----------------回代法求解L'x=y-----------------------------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('平方根法的解即为');function b=GJCholesky(A,b)n=size(A);n=n(1);v=zeros(n,1);%----------------------LDL'分解-----------------------------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);endB=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;end%-------------------前代法---------------------------A=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('改良平方根法解得的解即为');实验三、二次多项式拟合一、实验要求:用运算机语言编制利用QR分解求解线性最小二乘问题的通用子程序,用编写的程序求解一个二次多项式使在残向量的范数最小的意义下拟合下面的数据t-1 0iyi二、算法描述及实验步骤:QR分解求解程序如下:一、求A 的QR 分解; 二、计算b c 11T=Q ;3、求解上三角方程1c x =R 得x ;三、调试进程及实验结果:>> t=[-1 0 ]; >> y=[ ]; >> plot(t,y,'r*');>> legend('实验数据(ti,yi)'); >> xlabel('t'), ylabel('y');>> title('二次多项式拟合的数据点(ti,yi)的散点图');运行后屏幕显示数据的散点图(略).(3)编写以下MATLAB 程序计算)(x f 在),(i i y x 处的函数值,即输入程序 >> syms a b c >> t=[-1 0 ]; >> fi=a.*t.^2+ b.*t+c%运行后屏幕显示关于 ,,a b c 的线性方程组fi =[a-b+c,9/16*a-3/4*b+c,1/4*a-1/2*b+c,c,1/16*a+1/4*b+c,1/4*a+1/2*b+c,9/16*a+3/4*b +c]编写构造残向量2范数的MATLAB 程序>> y=[ ]; >> y=[ ];>> fy=fi-y; fy2=fy.^2; J=sum(fy.^2); 运行后屏幕显示误差平方和如下 J=(a-b+c-1)^2+(9/16*a-3/4*b+c-13/16)^2+(1/4*a-1/2*b+c-3/4)^2+(c-1)^2+(1/16*a+1/4*b+c-21/16)^2+(1/4*a+1/2*b+c-7/4)^2+(9/16*a+3/4*b+c-37/16)^2为求,,a b c 使J 达到最小,只需利用极值的必要条件0J a ∂=∂,0J b ∂=∂,0Jc∂=∂,取得关于,,a b c 的线性方程组,这能够由下面的MATLAB 程序完成,即输入程序 >> Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c);>> Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3) 运行后屏幕显示J 别离对,,a b c 的偏导数如下 Ja11 =451/128*a-63/32*b+43/8*c-887/128 Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*c-143/8解线性方程组112131000Ja Ja Ja ===,,,输入以下程序 >> A=[451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14]; >> B=[887/128,61/32,143/8]; >> C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f 及其系数C 如下 C =f =924/2999*x^2+10301/11996*x+4204/2999 故所求的拟合曲线为2()0.30810.8581 1.4018f x x x =++四、源程序:>> t=[-1 0 ]; >> y=[ ]; >> plot(t,y,'r*');>> legend('实验数据(ti,yi)'); >> xlabel('t'), ylabel('y');>> title('二次多项式拟合的数据点(ti,yi)的散点图'); >> syms a b c>> t=[-1 0 ]; >> fi=a.*t.^2+ b.*t+c fi =[ a-b+c, 9/16*a-3/4*b+c, 1/4*a-1/2*b+c, c, 1/16*a+1/4*b+c, 1/4*a+1/2*b+c, 9/16*a+3/4*b+c]>> y=[ ]; >> y=[ ];>> fy=fi-y; fy2=fy.^2; J=sum(fy.^2) J =(a-b+c-1)^2+(9/16*a-3/4*b+c-13/16)^2+(1/4*a-1/2*b+c-3/4)^2+(c-1)^2+(1/16*a+1/4*b+c-21/16)^2+(1/4*a+1/2*b+c-7/4)^2+(9/16*a+3/4*b+c-37/16)^2>> Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c);>> Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3)Ja11 =451/128*a-63/32*b+43/8*c-887/128Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*c-143/8>> A=[451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14]; >> B=[887/128,61/32,143/8];>> C=B/A, f=poly2sym(C)C =f =924/2999*x^2+10301/11996*x+4204/2999>>。

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

上机习题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代码%算法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));%算法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;end%高斯消去法解线性方程组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 ):由图,显然列主元消去法与精确解更为接近,不选主元的Gauss 消去法误差比列主元消去法大,且不如列主元消去法稳定。

Gauss 消去法重点在于A 的分解过程,无论A 如何分解,后面两步的运算过程不变。

2.先用你所熟悉的的计算机语言将平方根法和改进的平方根法编写成通用的子程序;然后用你编写的程序求解对称正定方程组Ax=b 。

Sol :(1)先用matlab 将平方根法和改进的平方根法编写成通用的子程序,得到L ,(D): 平方根法:L=Cholesky(A)改进的平方根法:[L,D]=LDLt(A)(2)求解得b Ly =求解得y x DL or y x L TT ==求解程序为x=Gauss(A,b,L,U,P)(TT DL U or L U == ,P 此时缺省,缺省时默认为单位矩阵)(3)计算脚本为ex1_2代码%算法function L=Cholesky(A)n=length(A);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:nA(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endendL=tril(A);end%计算LDL ‘分解:改进的平方根法function [L,D]=LDLt(A)n=length(A);for j=1:nfor i=1:nv(i,1)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1,1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1,1))/A(j,j);L=tril(A);D=diag(diag(A));L=L-diag(diag(L))+diag(ones(1,n));end%高斯消去法解线性方程组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;end%第二题%第一问A=10*eye(100)+diag(ones(1,99),-1)+diag(ones(1,99),1); b=round(100*rand(100,1));%平方根法L=Cholesky(A);x1_2_1_1=Gauss(A,b,L,L');%改进的平方根法[L,D]=LDLt(A);x1_2_1_2=Gauss(A,b,L,D*L');%第二问A=hilb(40);b=sum(A);b=b';%平方根法L=Cholesky(A);x1_2_2_1=Gauss(A,b,L,L');%改进的平方根法[L,D]=LDLt(A);x1_2_2_2=Gauss(A,b,L,D*L');结果分别为x1_2_1_1 =7.25868.4143-0.40138.59840.2249 2.3336 4.4389 8.2772 8.7890 -0.1667 8.8784 8.3824 3.2978 7.6401 0.3014 3.3457 8.2418 6.2368 8.3906 5.8575 -0.9656 7.7981 7.98425.36016.4143 6.4966 2.6201 6.3024 0.35637.13422.8506 0.1927 0.2227 7.5801 5.9762 1.6583 9.4409 -1.0677 4.2362 2.70606.70377.2570 0.7265 4.4778 3.4959 5.56375.86756.7614 1.5180 6.0582 5.9002 0.9397 0.7031 4.0294 9.0032 1.9382文档来源为:从网络收集整理.word版本可编辑.欢迎下载支持.5.61500.91207.26521.43604.37495.81467.47918.39424.57890.81691.25231.66038.14480.89157.94010.70758.98492.44371.57771.77905.63193.90182.35067.59254.72454.16278.64836.8087 6.5589 2.6027 5.4140 0.2577 0.0090 4.6522 6.4685 8.6626 -0.0948 5.2856 4.2385 -0.67063.4671 x1_2_1_2 =7.25868.4143 -0.4013 8.5984 5.4177 0.2249 2.3336 4.4389 8.2772 8.7890 -0.16678.3824 3.2978 7.6401 0.3014 3.3457 8.2418 6.2368 8.3906 5.8575 -0.9656 7.7981 7.98425.36016.4143 6.4966 2.6201 6.3024 0.35637.1342 -0.6985 2.8506 0.1927 0.2227 7.5801 5.9762 1.6583-1.0677 4.2362 2.70606.70377.2570 0.7265 4.4778 3.4959 5.56375.86756.7614 1.5180 6.0582 5.9002 0.9397 0.7031 4.0294 9.0032 1.9382 5.6150 0.91207.2652 1.43604.37495.8146 7.47914.57890.81691.2523 1.6603 8.1448 0.8915 7.9401 0.7075 8.98492.4437 1.5777 1.7790 5.63193.9018 2.3506 7.59254.7245 4.1627 8.6483 1.3543 6.8087 6.5589 2.60275.4140 0.2577 0.00906.4685 8.6626 -0.0948 5.2856 4.2385 -0.67063.4671x1_2_2_1 = 1.0e+07 * 0.0000 -0.0000 0.0001 -0.0004 -0.0014 0.0424 -0.2980 1.1419 -2.7335 4.2539 -4.3018 2.7733 -1.1989 0.5406 -0.3688 0.3285 -0.4438-0.2513 0.0565 0.0000 -0.0051 0.0071 -0.0027 -0.0031 0.0036 -0.0019 0.0009 0.0002 -0.0002 -0.0006 0.0004 0.0001 -0.0002 0.0001 0.0000 -0.0000 0.0000 -0.0000 -0.0000 x1_2_2_2 = 1.0000 1.0000 0.99981.00640.86811.8034 -1.5693 5.5763 -2.5315 -1.7693 10.4883 -6.2807 0.5882 -4.7157 22.8299 -19.9134 8.7032 10.3265 -25.2140 10.0282 12.3882 -1.9425 14.1891 -12.0552 -0.5803 -12.47918.56529.8724 -10.5502-5.813213.421611.1767-64.315446.383712.6957-21.755612.1204-1.93423.用第1题的程序求解第2题的两个方程组并比较所有的计算结果,然后评价各个方法的优劣。

相关文档
最新文档