五对角稀疏矩阵求解.docVIP

第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=2 U(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-2 i=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=n U(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分解的比较 N LU(cpu time ) sparselu(cpu time) 10 2.7800e-05 0.000956 100 0.000105 0.006699 1000 0.001014 0.072862 10000