五对角稀疏矩阵求解.doc

第01次大作业课程名称： 计算方法 实验项目名称： Solve a five-diagonal linear system(Topic 1)学生姓名： 张弘 一、题目二、问题分析及背景1)稀疏矩阵求解的应用：矩阵求解是数值计算的核心，而稀疏矩阵的求解则是数值计算的关键步骤。而时间及空间（内存）的不足往往成为矩阵求解的瓶颈。对于方程组求解一般有一下三种方式。1.迭代法：迭代法是通过给定初始向量x0，用设计好的迭代公式计算出向量序列x1, x2,x3 xk，直至收敛到方程组的解。其特点是不改变矩阵非零元的位置，可以保持牺牲矩阵的稀疏性，适合求解大型稀疏线性方程组。常用的迭代法有：Jacobi迭代法，G-S迭代法，SOR迭代法等2.直解法：直接法是在不考虑误差的前提下，通过有限步的四则运算求的方程组解的方法。其特点是计算精度高，但计算需要较多存储空间，计算量大。其适合于求解规模不大的线性方程组及非零元集中在对角元附近的大型稀疏方程组。常用的方法有：LU分解(高斯消去法)，Cholesky分解法,frontal， multifrontal(多波前法)， supermodel分解等。其中LU分解法相比其他解法有以下优势：保持稀疏性，较高稳定性，可充分发挥cpu效率。其不足是算法复杂。对于稀疏的对角方程组可以使用追赶法进行LU分解，这种方法可以大量减少运算提高运算效率。3.多重网格法多重网格方法是由对偏微分方程里得出的代数方程组的求解的研究引发出来的一种计算方法，它已经成为求解大型科学与工程计算问题的最有效方法之一。多重网格方法包括两种不同类型方法，根据构造方法（空间和基函数）的不同,分为一般多重网格方法和多水平方法。一些论著中也将这两种方法统称为多水平方法三、算法设计及分析1）因位我们是针对五对角矩阵进行求解，由于这样的矩阵是稀疏的，观察题目中给出的例子也可以看出其分解后的矩阵L和U都是稀疏矩阵，因此参考教材第21页的lin_band2.m程序，按如下方法对五对角矩阵A做LU分解。首先A的LU分解有如下形式：因此利用比较系数法得到如下递推形式：对于i=1：U(1,1) = A(1,1);L(2,1) = A(2,1)/U(1,1);L(3,1) = A(3,1)/U(1,1);I=2U(1,2) = A(1,2);U(2,2) = A(2,2) - L(2,1)*U(1,2);L(3,2) = (A(3,2) - L(3,1)*U(1,2)/U(2,2);L(4,2) = A(4,2)/U(2,2); U(i-1,i) = A(i-1,i) + L(i-1,i-2); U(i,i) = A(i,i) - L(i,i-2)*A(i-2,i) - L(i,i-1)*U(i-1,i); L(i+1,i) = (A(i+1,i) - L(i+1,i-1)*U(i-1,i)/U(i,i);L(i+2,i) = A(i+2,i)/U(i,i);i=3 4 .n-2i=n-1时U(n-2,n-1) = A(n-2,n-1) + L(n-2,n-3);U(n-1,n-1) = A(n-1,n-1) + L(n-1,n-3) - L(n-1,n-2)*U(n-2,n-1);L(n,n-1) = (A(n,n-1) - L(n,n-2)*U(n-2,n-1)/U(n-1,n-1);I=nU(n-1,n) = A(n-1,n) + L(n-1,n-2);U(n,n) = A(n,n) - L(n,n-2)*A(n-2,n)-L(n,n-1)*U(n-1,n);因此可根据以上递推形式写出程序。由于LU矩阵是稀疏的，所以可以做如下处理：令L U均为N*2维的向量，这样就可以只存储需要我们计算的L 和U中的元素，而不必将其以矩阵形式存储。四、实验结果比较1）通过利用matlab程序自带的标准LU分解函数lu()，稀疏矩阵的LU分解函数ilu(), 和按照以上算法实现的sparselu()，通过对n=10 100 1000 10000 100000 100000 阶矩阵进行LU分解比较，可以得出如下结果：1.matlab的lu()函数与sparselu()分解的比较。表1.matlabLU分解和 sparselu()进行LU分解的比较NLU(cpu time )sparselu(cpu time)102.7800e-050.0009561000.0001050.006699 10000.0010140.072862 100000.0099681.1213151000000.10998347.0531810000001.101775inf 分析1： 由上表比较可知，随着矩阵规模N每增加一个数量级 ，matlab的lu函数及自己的sparselu运行时间均相应的增加一个数量级。且可以看出，虽然利用追赶法求解对5对角系数矩阵进行LU分解可以减少运算量，当由于matlab内部的并行机制相当好，分解相同规模的五对角稀疏矩阵matlab lu函数的运行时间大致比sparselu的运行时间快一至两个个数量级。 2)稀疏矩阵直接lu分解法求五对角方程组，matlab左除法”和GaussSeidel迭代法求解的精度与cpu时间比较。表2.matlab sparsesolve() 与GaussSeidel迭代法求解方程的比较NSparsesolveTimeErrorTimeError100.0000131.3911e-150.0006281.1854e-151000.0000583.2448e-130.0072864.4236e-1310000.0003731.1354e-100.0799331.0811e-10100000.0036643.5262e-081.1667893.4672e-081000000.0434731.0578e-0547.9028181.0430e-0510000000.4418130.0034-表示时间太长，超过10min表3.要求不同精度时，GaussSeidel迭代法所需要的时间NTime1/sError1Time2/sError2100.00409.1399e-110.05607.6113e-161000.27589.9993e-064.36463.4632e-13100019.961711.2279-10000-100000-1000000-当5对角方程组的规模大于1000时，虽然GaussSeidel迭代法依然是收敛的，但由于其收敛速度相当慢，最初迭代时收敛稍微快一点，担当迭代精度趋近10-102时每迭代一次精度下降很少几乎不足0.01，而且规模愈大精度下降愈慢.在规模为1000时学迭代10000次精度才达到11.2279，因此未对10000 1000000规模下的方程迭代求解，上表中也未列出方程组规模为10000 1000000 的情况下其收敛情况。 由上表比较可知，在计算效率上来讲matlab的左除法效率最高，能求解的矩阵规模最大，其次是sparsesolve可以求解105阶矩阵，GaussSeidel迭代法最慢，只可以求解103阶矩阵。并且有GaussSeidel迭代法运行时间比Sparsesolve高一个数量级，sparsesolve比高一个数量级。在精度上由比较知：虽然GaussSeidel，sparsesolve与matlab的相比计算效率低但三中方法计算精度相当。 五、实验总结这次试验通过采取直接LU分解对不同规模的五对角方程组进行求解，并将求解结果与matlab自带的函数lu(), ilu(), 以及GaussSeidel迭代法求解方程进行了比较。最后得到这样的结论：matlab的函数lu() 及除法 针对稀疏矩阵进行了优化，使其可以在求解大规模及超大规模的稀疏矩阵方程组时能有很好的效果。其次，在试验中还可以发现，针对此次作业中的五对角方程组而设计的直接LU分解（追赶法）虽然相对于标准LU分解的运算少了很多，但由于程序优化不如matlab的lu分解，追赶法的运行效率相对较低其精度与matlab自带的左除法求解精度相当的，且要比GaussSeidel迭代法效率高。因此具体问题具体分析有时候可以提高求解的质量，当然这要牺牲一定的时间为代价。六、参考文献1叶兴德,程晓良,陈明飞,薛莲数值分析基础.浙江大学出版社附录：程序代码（sparselu.m comparelu.m gs_mt.m）1）comparesolu.m%Matlab的LU分解和sparselu(追赶法)的运行时间比较function compareluN=10 100 1000 10000 100000 1000000;%生成五对角稀疏矩阵%生成右端向量b = rand(n,1);L = zeros(n,2);U = zeros(n,2);%m = 100;%存储LU分解的时间t1=zeros(6,1);%存储sparselu分解的时间t2=zeros(6,1);for i =6 n = N(i); e = ones (n,1); A=spdiags(-e,-e,4*e,-e,-e,-2,-1,0,1,2,n,n); % lu分解 tic for j = 1:m L1,U1 = lu(A); end t1(i) = toc; % sparselu分解 for j=1:m L2,U2 = sparselu(A); end t2(i) = toc; endt1 = t1/m;t2 = t2/m;fprintf( N LU(cpu time ) sparselu(cpu time)n);for i=6fprintf(%6d %10.6f %10.6f n,N(i),t1(i),t2(i);end2)compareSolve.m%Matlab的和sparsesolve(追赶法)求解方程组的运行时间比较function compareSolveN=10 100 1000 10000 100000 1000000;%生成五对角稀疏矩阵%生成右端向量%m = 2;%存储LU分解的时间t1=zeros(6,1);%存储sparselu分解的时间t2=zeros(6,1);t3=zeros(6,1);err1=zeros(6,1);err2=zeros(6,1);err3=zeros(6,1);for i =3 n = N(i); e = ones (n,1); A=spdiags(-e,-e,4*e,-e,-e,-2,-1,0,1,2,n,n); b = rand(n,1); tic for j=1:m x2 = sparsesolve(A,b); err2(i) = err2(i) + norm(A*x2-b); end t2(i) = toc; err2(i) = err2(i)/m; tic for j=1:m x1 = Ab; err1(i) = err1(i) + norm(A*x1-b); end t1(i) = toc; err1(i) = err1(i)/m; tic for j= 1:m x3 = gs_mt(A,b,1e-1); err3(i) = err3(i) + norm(A*x3-b); end t3(i) = toc; err3(i) =err3(i)/m;endt3 =t3/m;t1 = t1/m;t2 = t2/m;fprintf( N matlab (cpu time )sparselu(cpu time) GaussSeideln);for i=1:6fprintf(%6d %10.6f %10.6f %10.6f n,N(i),t1(i),t2(i),t3(i);endfprintf(N Norm(A*x-b)matlab sparselu GaussSeideln);for i=1:6fprintf(%6d %10.6f %10.6f %10.6f n,N(i),err1(i),err2(i),err3(i);end3)sparselu.m%使用追赶法求解5对角系数矩阵function Lm,Um = sparselu(A)%5对角稀疏矩阵LU分解%调用方式%sparse(A) A为所求矩阵%Lm 和Um分别返回分解矩阵L U%矩阵A举例：%e = ones (n ,1);%A=spdiags(-e,-e,4*e,-e,-e,-2,-1,0,1,2,n,n);n = size(A,1);%以二维向量形式存储LUL = zeros(n,2);U = zeros(n,2);%开始对A做LU分解U(1,2) = A(1,1);U(2,1) = A(1,2);L(1,1) = A(2,1)/U(1,2);L(1,2) = A(3,1)/U(1,2);U(2,2) = A(2,2) - L(1,1)*U(2,1);L(2,1) = (A(3,2) - L(1,2)*U(2,1)/U(2,2);L(2,2) = A(4,2)/U(2,2);for i = 3:n-2 U(i,1) = A(i-1,i) + L(i-2,1); U(i,2) = A(i,i) - L(i-2,2)*A(i-2,i) - L(i-1,1)*U(i,1); L(i,1) = (A(i+1,i) - L(i-1,2)*U(i,1)/U(i,2); L(i,2) = A(i+2,i)/U(i,2);endU(n-1,1) = A(n-2,n-1) + L(n-3,1);U(n-1,2) = A(n-1,n-1) + L(n-3,2) - L(n-2,1)*U(n-1,1);L(n-1,1) = (A(n,n-1) - L(n-2,2)*U(n-1,1)/U(n-1,2);U(n,1) = A(n-1,n) + L(n-2,1);U(n,2) = A(n,n) - L(n-2,2)*A(n-2,n)-L(n-1,1)*U(n,1);%LU分解结束%将L U向量稀疏化为LU矩阵Lm和Ume = zeros(n,1);Lm = spdiags(L(1:n,2),L(1:n,1),e,-2,-1,0,n,n);Um = spdiags(U(1:n,2),U(2:n,1); 0,-e,0 1 2,n,n);4)sparsesolvefunction x = sparsesolve(A,b)%5对角稀疏矩阵追赶法法求解%调用方式 x=sparsesolve(A,b) A为系数矩阵 b为右端向量%矩阵A举例：%A=spdiags(-e,-e,4*e,-e,-e,-2,-1,0,1,2,n,n);n= size(A,1);e = ones (n ,1);%以二维向量形式存储LUL = zeros(n,2);U = zeros(n,2);%开始对A做LU分解U(1,2) = A(1,1);U(2,1) = A(1,2);L(1,1) = A(2,1)/U(1,2);L(1,2) = A(3,1)/U(1,2);U(2,2) = A(2,2) - L(1,1)*U(2,1);L(2,1) = (A(3,2) - L(1,2)*U(2,1)/U(2,2);L(2,2) = A(4,2)/U(2,2);for i = 3:n-2 U(i,1) = A(i-1,i) + L(i-2,1); U(i,2) = A(i,i) - L(i-2,2)*A(i-2,i) - L(i-1,1)*U(i,1); L(i,1) = (A(i+1,i) - L(i-1,2)*U(i,1)/U(i,2); L(i,2) = A(i+2,i)/U(i,2);endU(n-1,1) = A(n-2,n-1) + L(n-3,1);U(n-