安徽工业大学数理科学与工程学院
病态线性方程组的求解
专业数学与应用数学
班级数***班
学号*******
姓名 ***
指导教师***
二○一五年五月
一、设计目的:
为了更加透彻的熟悉数值分析课程,学习各种数学软件的使用,锻炼自己对知识的实际运用能力。
二、引言:
用直接法求解AX=F 线性方程组,对于系数矩阵A 对角占优是很有效的。
方程阶数不高时,人们经常使用;而当方程组阶数大时,由于积累误差,导致结果失真。
为了克服误差积累问题,通常用迭代法。
它具有可达到所要求的精度和对计算内存要求不大的优点,对求解大型线性方程组,迭代法计算时间远比直接法少,所以在实际计算中,迭代法也被人们广泛使用。
然而迭代法要研究迭代格式的收敛性,如Jacobi 迭代对系数矩阵为病态矩阵不收敛,为此我们提供一种修改的Jacobi 迭代,并给出一些数值例子来说明有较好的效果。
三、解线性方程组迭代法的描述
设线性方程组AX=P ,这里A:{a ij },X:{x i },F:{f i },1≤i,j ≤n,为了更广泛地应用,对A 只限制实的非奇异矩阵,那么,若给定初值x )0(,我们熟知的有: Jacobi 迭代:
x )
1(+k i
=(f i -∑n
j
k j ij x a )()/a ii j i ≠, 1≤i ≤n
四、求解病态线性方程组的另一种迭代解法
设线性方程BX=F, 这里系数矩阵B 是病态的,指的是矩阵条件数是较大的。
条件数越大,就越难求得准确解,为此,我们将方程的两端同加DX 项,那 么相应的Jacobi 迭代有:
X ])([)()(1)1(k k X H D F A D -++=-+ (1)
这里,A 为B 的对角阵,即A: {b ii },H:{b ij } j i ≠ 1≤i,j ≤n,
记M= )()(1H D A D -+-,)()1()(k k k X X X -=∆+,那么有:
)0()1()(X M X M X k k k ∆==∆=∆-
由此可见(1)式收敛,迭代矩阵M 的谱半径应满足p(M )<1,谱半径若用M 的特征值判断,则较为繁锁;若B 不可约,根据线性代数对角占优简单迭代必收敛的性质,为简单起见,我们取D 为对角线阵,即D: {d i },为保证收敛就得取d i =Sign{∑j
ii ij b b |,
|},符号Sign{a,b}的含义是与b 同号,数值
取a,这是充分条件,实际计算中有时可放宽处理,比如可取d i =Sign{ii ij ij
b b Max |,
|}或者d i =Sign{ii ij j
b b Max |,
|} j i ≠,因而相应
的Jacobi 迭代修改为:
x )1(+k i =(f i -∑≠+i
j k i i k j ij x d x b )()()/(b ii +d i ) (2)
下面列举普通Jacobi 迭代不收敛,解不出正确结果,而用修改的Jacobi 迭代可求出满意结果的例子:
例 1:
A=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--211111112 F=⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡--032 精确解为X=[]T
111---,普通Jacobi 迭代不收敛,取
d i =Sign{ii ij j
a a Max |,
|} j i ≠, 1≤i,j ≤3
用(2)式迭代40步,X=[]T 000003.1000053.1000035.1---
例 2:
方程系数为Hilbert 阵
H:{)1/(1-+=j i h ij } F:{∑=j
ij i h f } 1≤i,j ≤n
精确解X=[]T 111 ,普通Jacobi 迭代发散,取d i =Sign{∑j
ii ij h h |,
|}用
修改的(2)式迭代,n=1200时迭代10870步结果摘录如下:
五、结论与问题
以上数值试验表明该迭代算法对求解病态线性方程组是有效的,其优点是可达到预定的精度。
求解病态方程组大条件数的系数矩阵,要获得较正确的解是很困难的。
本文的算法迭代虽能保证收敛,但与精确解的误差到底怎样,还应将解代入原方程,衡量、检验求得解的可信程度;另外,如何克服求解病态线性方程组收敛缓慢的问题,还有待于进一步工作。
附录
例 1的matlab计算程序:
clear
m=40;
n=3;
A=[2 -1 1;1 1 1;1 1 -2];
F=[-2 -3 0];
x=zeros(n,m);
for i=1:n
x(i,1)=0;
end
for i=1:n
if i==1
d(i)=max(abs(A(i,2:n)))*A(i,i)/abs(A(i,i));
elseif i==n
d(i)=max(abs(A(i,1:n-1)))*A(i,i)/abs(A(i,i));
else
d(i)=max(max(abs(A(i,1:i-1))),max(abs(A(i,i-1:n))))*A(i
,i)/abs(A(i,i));
end
end
for k=2:m
for i=1:n
x(i,k)=(F(i)-A(i,:)*x(:,k-1)+A(i,i)*x(i,k-1)+d(i)*x(i,k-1))/ (A(i,i)+d(i));
end
end
x(:,m)
例 2的matlab计算程序:
clear
n=1200;
m=10870;
for i=1:n
for j=1:n
H(i,j)=1/(i+j-1);
end
end
for i=1:n
f(i)=sum(H(i,:));
end
for i=1:n
d(i)=sum(abs(H(i,:)))*H(i,i)/abs(H(i,i));
end
for i=1:n
x(i,1)=0;
end
for k=1:m
for i=1:n
x(i,k+1)=(f(i)-H(i,:)*x(:,k)+H(i,i)*x(i,k)+d(i)*x(i,k))/( H(i,i)+d(i));
end
end
x(:,m)。