数值分析课件典型例题与习题2

1、三、四章内容提要 典型例题分析,数值分析典型例题 II, ,11:38,化难为易,化繁为简,化繁为简,初等行变换不改变方程组的解,1.交换矩阵的第i行与第j行的位置,2.以非零数k乘以矩阵的第i行的每个元素,3.把矩阵的第i行的每个元素的k倍加到第j行的对应元素上去,A(n 1) = Fn-1Fn-2F1 A,其中Fk 为 Frobenius矩阵。,A=F1-1F2-1 Fn-1-1 A(n 1),直接方法: 高斯消元法,L U,高斯消元法本质是矩阵的分解。,矩阵分解(Top 10 Algorithms),(1)特征值分解: A=CDC, C,D=eig(A),(3) LU分解: PA=LU,

2、 L,U,P=lu(A),(4)Cholesky分解: A=LLT, R=chol(A),11:38,(2)奇异值分解: A=USV , U,S,V = svd(A),(5) 非负矩阵分解 Learning the Parts of Objects by Non-negative Matrix Factorization, Nature, 1999,Demo1,I=imread(monalisa.pgm); U,S,V=svd(double(I); s=diag(S); n1=5; Snew=diag(s(1:n1);zeros(size(s,1)-n1,1); figure,imshow(U

3、*Snew*V,) n2=20; Snew=diag(s(1:n2);zeros(size(s,1)-n2,1); figure, imshow(U*Snew*V,),11:38,迭代法思想:,11:38,Iterate:To say or doagainoragain and again,迭代背后的思想是一种与传统思维模式截然不同的方式,传统思维方式往往希望一遍做好,一次成功;但是迭代开发意味着反复地做,不断地根据反馈进行调整。,11:38,迭代格式构造,收敛条件(局部vs全局),中止准则,加速(松弛思想),Aitken加速方法 超松弛加速方法,11:38,共轭梯度法的关键是构造一组两两共轭

4、的方向(第k步迭代生成共轭方向张成k维子空间)。巧妙的是共轭方向可以由上次搜索方向和当前的梯度方向组合产生。,现代迭代方法 (Top 10 Algorithms),最速下降法思想简单,但收敛速度慢。本质上因为负梯度方向函数下降快是局部性质。,复杂性:,高斯消元法共用乘法和除法次数为n3/3+ n2-n/3,常用记号O表示是多少阶的,则O (n3/3)。,注释2: 复杂性对估计求解大型方程组所需的时间有用。例如在一台特定的计算机上求解n=500个方程的方程组所需的时间我们可以通过求解一个n=50个方程的方程组得到一个很好的猜测,即对用掉的时间按比例放大1000倍。,注释1:按O的记法,把n的不同

5、次幂相加的结果仅保留了最高次幂,因为最高次幂决定了当n趋近无穷时的极限形态。换而言之,对于大的n ,低阶项对算法的执行时间的估计没有太大影响, 仅需要近似估计执行时间时可以忽略不计。,常用的范数:,Hilbert矩阵条件数: for i=1:10 c(i)=cond(hilb(i),2);%vander(1:i) end,plot(1:10,c),范数(全局),问题的好与坏,算法的快与慢,| X(k+1) X* | |B| | X(k) X* |,范数的威力和魅力:,ekT= 0 0 1 0 0 ,= I mkekT,Fk1 = I + mk ekT,( k = 1, 2, , n 1),F1

6、-1F2-1 Fn-1-1 =, ,Fk1 Fj1 = I + mk ekT+mj ejT kj,F1-1F2-1 Fn-1-1 = I + m1 e1T+mn-1 en-1T,例1.,ekTmj =0, kj,例2.设A为对称矩阵。高斯消元法一步后,A约化为,证明 B 也是对称矩阵。,约化主元不为零的判断,定理3.1 约化主元 的充分必要条件是矩阵A各阶顺序主子式 不等于零。即,例4.严格对角占优矩阵的约化主元ak,k(k-1) 0 (k=1,n) 。,例3.严格主对角占优矩阵一定是非奇异的。,严格对角占优矩阵:高斯消元法中约化主元不等于零,Jacobi方法, GS方法和SOR方法收敛。,思

7、考: 若A是严格对角占优矩阵, 当0 w =1时, SOR方法收敛。,例5.Ax=b,其中A对称正定,问解此方程的Jacobi迭代法是否一定收敛？,对称正定矩阵:直接法高斯消元法中约化主元大于零,迭代法GS方法,SOR方法,最速下降方法和共轭梯度法收敛。,例6.,例7.,例8.,病因是条件数大,病症是什么呢?,例9.矩阵的Doolittle分解,A = LU, L为单位下三角矩阵,U为上三角矩阵。,a11= u11, , a1n= u1n,a21 = m21u11, , an1 = mn1u11,11:38,对A的元素aij ,当 jk 和 ik+1时,m21u12+ u22=a22, , m

8、21u1n+ u2n=a2n,u22=a22 - m21u12, , u2n=a2n- m21u1n,m31u12+ m32u22=a32, , mn1u12+ mn2u22=an2,m32=(a32- m31u12)/u22, , mn2=(an2- mn1u12)/u22,11:38,矩阵L和矩阵U的元素计算公式,计算次序,11:38,更新顺序: 先行后列 列除行不除 旧元素减去所在行和列前k-1分量点乘的加和,Doolittle分解：,更新顺序: 先列后行 行除列不除 旧元素减去所在行和列前k-1分量点乘的加和,Crout分解：,求矩阵的Doolittle分解,11:38,三对角矩阵分解

9、,11:38,步骤I:,三对角矩阵分解 A = L U,( k = 2, 3, , n ),11:38,下三角方程组 LY = f,步骤II:,上三角方程组 UX = Y,步骤III:,11:38,求解三对角方程组的追赶法,11:38,对称正定矩阵的 Cholesky分解,(1) 对于非零向量, xTAx总是正的; (2) A的顺序主子式全大于零; (3) A的特征值全为正实数。 help chol,对称正定矩阵:,1) 序列收敛 2) 迭代矩阵谱半径小于1 3) 迭代矩阵特征值全小于1 4) 迭代矩阵范数小于1 5) 反证法,迭代法收敛证明思路:,例10.设A是一个可逆矩阵, 矩阵序列满足

10、Xk+1=Xk(2I A Xk )，（k =0,1,2,） 则当 时,证明：由Xk+1=Xk(2I A Xk )，得 I AXk+1 = I A Xk(2I A Xk )= (I A Xk )2 于是 I AXk =(I A Xk -1)2 =(I A Xk -2)22 = ,例11.,思路:,(1) A1 = B ( I + R + R2 + )； (2)任意给定n阶矩阵X0,由迭代格式 Xk+1 = Xk R + B ( k = 0,1,2, ) 产生的矩阵序列 Xk 收敛到矩阵A-1； (3)对矩阵序列 Xk , 有误差估计式,例12.设A是n阶可逆矩阵,有A的一个近似逆B,令 R=I

11、AB。如果 | R | q 1,试证明,收敛的w取值范围。,例13. 方程组Ax=b, 其中A是对称正定阵,讨论迭代格式,思路:,例14.,思路:,续例14.,思路:,例15.,思路:,定理4.1 对任意的f和任意的初始向量X(0)迭代法 X(k+1) =B X(k) +f 收敛的充分必要条件是,例16.,思路:,例17.,思路: -1/2a1/2收敛, -1/2a1时系数矩阵正定。,例18.,例18.,例18.,参考:Writting Fast Matlab Codes,1.中小规模线性方程组 x = Ab; % Solves A*x = b,2.(超)大规模线性方程组(Preconditi

12、oned cg),作业,题目1: 从理论角度(复杂度和收敛性)比较各 种方法的优劣。,11:38,题目2: 从数值角度比较各种方法的性能(公平),题目3: 科研或生活中遇到的线性方程组？,提示:,参考代码: Matlab代码Iterative Solver,11:38,1. 如何生成方程组 A = gallery(poisson, n); b=ones(size(A,1),1); x0=zeros(size(A,1),1);,2. 如何生成方程组 矩阵市场: /MatrixMarket/ UF Sparse Matrix Collection: http:

13、//research/sparse/matrices/,提示:,11:38,3. 更具挑战和趣味的例子 图像编辑 参考: Matlab代码Possion 图像复原 参考: Deblurring Images: Matrices, Spectra, and Filtering Google PageRank 参考: Numerical Computing with MATLAB 数据挖掘和模式识别 参考: Matrix Methods in Data Mining and Pattern Recognition 高光谱图像解混 参考: Sparse Unmixing of Hyperspect