稀疏矩阵在迭代求解线性方程组中的运用

稀疏矩阵在迭代求解线性方程组中的运用

学院：自动化院

专业：电力系统及其自动化

姓名：张庆磊

学号：111101112

指导老师：杨伟

摘要:对于稀疏矩阵的稀疏存储技术进行了研究，研究了按行存储的检索方式，以便于迭代计算，分别对雅克比法和高斯—赛特法的算法进行理论研究和程序实现，并且比较了两种方法的优劣。 关键词：系数技术；线性方程；迭代

大型稀疏矩阵线性化方程组的求解问题，在电力系统中有着广泛的的运用。由于电力网本身的结构限制，节点导纳矩阵节点繁多，而仅有少量的非零元，稀疏度很高，若采用传统存储计算方式，会占用大量的存储空间，并且降低运算效率。在迭代计算中，由于无法分辨零元素，也会无谓地浪费运算时间。因此稀疏技术在求解方程组中的运用显得尤为重要。

1.稀疏矢量与稀疏矩阵的存储

稀疏矢量与稀疏矩阵的存储特点是排零存储,即只存储其中的非零元和有关的检索信息.存储的目的是为了在计算中方便的访问和运用，这就要求既节省内存，又便于搜索。论文采用了按行存储格式。

按行顺序依次存储A 中的非零元，同一行元素依次排列在一起，存储格式： V A ——按行存储矩阵A 中的非零元ij a ，共τ个； JA ——按行存储矩阵A 中的非零元的列号，共τ个； IA ——记录A 中每行第一个非零元在V A 中的位置，共n 个。

2.迭代法求解线性方程

设n n

n n

R b R A ∈∈⨯,,考虑线性方程组 Ax=b

一般的，先将式1变为同解方程组 X=Bx+f 形成迭代式

f Bx x k k +=+)()1(

式中：B 为迭代矩阵；

若要求迭代式收敛则需满足1B 1)(<<或B ρ，其中∙为任意范数 式1中，A 可以分裂为 A=M+N

(1)

其中，M 非奇异，则可以得到

b M Nx M x 11--+-=

(2)

令

A M I N M

B 11---=-= (3) b M f 1-=

(4)

对于以矩阵A ，有A=D+L+U

其中D 为对角线矩阵，L 为严格下三角矩阵，U 为严格上三角矩阵 由此构造迭代法，令M=D ，N=L+U

f Bx x k k +=+)()1(

(5)

式中，向量f 和迭代矩阵B 为

⎪⎩⎪⎨⎧-=+-==---A

D I U L D B b

D f 1

11

)( (6)

上式称为Jacobi 迭代。 可简单描述为

ii

n

i j j k j ij i k i

a x a

b x

1,1)()1(*⎪⎪⎭⎫ ⎝⎛-=∑≠=+ (7)

如果令M=D+L ，N=U ，对应的分裂式有

⎪⎩⎪⎨⎧+-=+-=+=---A

L D I U L D B b

L D f 1

11

)()()( (8)

便得到Gauss —Seidel 迭代法，即

ii

n

i j k j

ij

i j k j ij i k i

a x

a

x a b x

11

)

(11)

1()1(*⎪⎪⎭⎫

⎝

⎛--=∑∑+=-=++ (9)

3.J 迭代法和G —S 迭代法的收敛性

若n n ij R a A ⨯∈=)(满足

),,2,1(,1n i a

a n

i

j j ij

ii ∑≠==>

(10)

则称A 为严格对角占优矩阵。

若A 为严格对角占优矩阵，则Ax=b 的J 迭代法和G —S 迭代法均收敛。但特别的是，J 迭代法和G —S 迭代法的收敛没有相容性，即J 迭代法下收敛的矩阵在G —S 迭代法上无法确定是否收敛，反之亦然。

4.算法程序的设计

流程图：

开始

输入A 、b 、o

将A 使用稀疏技术存储

O=1

Jacobi 迭代

Gauss-Seidel

迭代

误差>e

误差>e

输出结果以及迭代次数

结束

Y

N

Y

Y

N

N

具体程序：

A=input('Enter Matrix A=')

b=input('Enter Matrix B=')

b=input('Enter n=')

int i;

int j;

int ia;

int ja;

int vl;

int vu;

int vk;

int vd; %%vl:存储下三角 vu:存储上三角 vd:对角元素 vk:对应行号

[ia,ja]=size(A);

VA=zeros(3,ja);

vu=1;

for i=1:ia

vk=0;

for j=1:ja

if((A(i,j)~=0)&(vk==0))

VA(3,i)=vu;

end

if(A(i,j)~=0)

vk=1;

VA(1,vu)=A(i,j);

VA(2,vu)=j;

vu=vu+1;

end

end

end

VA(3,ia+1)=vu;

VA(3,ia+2)=vu;

vs=1;

X=ones(1,ja)

s=2

while((s>0.0001)&(vs<100))

vs=vs+1

s=0

for i=1:ia

for j=VA(3,i):(VA(3,i+1)-1)

if(VA(2,j)==i)

vd=j

end

end