数值线性代数课设课件资料.doc

数值线性代数课程设计报告

姓名:陶英
学号:081410124

任课教师:杨熙

南京航空航天大学
2016 年 6 月 22日

求解线性方程组的三种迭代法及其结果比较

摘要
当今的环境下，数值计算越来越依赖于计算机。大规模科学计算和工程技术中许多问题的解决，最终归结为大型稀疏线性方程组的求解，其求解时间在整个问题求解时间中占有很大的比重，有的甚至达到80%。由于现今科学研究和大型项目中各种复杂的可以对计算精度和计算速度的要求越来越高。因此，作为大规模科学计算基础的线性代数方程组的高效数值求解引起了人们的普遍关注。这种方程组的求解一般采用迭代法。
关于迭代法，是有很多种解决公式的：Jacobi，G-S和超松弛迭代法。这三种方法的原理大致相同，Jacobi需要给定初向量，G-S则需要给定初值，超松弛法是对Guass-Seidel迭代法的加权平均改造。而本文则是对大型稀疏线性方程组迭代求解与三种迭代法（Jacobi，Gauss-Seidel和超松弛迭代法）的收敛速度与精确解的误差比较做出研究。

关键词：Jacobi迭代法；Gauss-Seidel迭代法；SOR迭代法；线性方程组

1 方法与理论的叙述
1.1迭代法简介
1.Jacobi迭代法：
对于非奇异线性方程组Ax=b，令A=D-L-U，其中

则原方程组可改写为：
（2.2）
其中

给定初始向量：

由（2.2）可以构造迭代公式：

其分量形式为：

2. Guass-Seidel迭代法：
类似于Jacobi迭代法，给定初值：

令

则得到Guass-Seidel公式：

其分量形式为：

3.超松弛迭代法（SOR 迭代法）：
SOR迭代法是对Guass-Seidel迭代法的加权平均改造，即

为Guass-Seidel迭代解，即

它的分量形式为：

其中ω称为松弛因子，当ω>1时称为超松弛；当ω<1时叫低松弛；ω=1时就是
Guass-Seidel迭代。

上述三种经典迭代法收敛的充分必要条件是迭代矩阵谱半径小于1。
谱半径不易求解，而在一定条件下，通过系数矩阵A的性质可判断迭代法的收敛性。
定理1：
若系数矩阵A是严格对角占优或不可约对角占优，则Jacobi迭代法和Gauss-Seidel迭代法均收敛。
定理2：
（1）SOR迭代法收敛的必要条件是0<w<2；
（2）若系数矩阵A严格对角占优或不可约对角占优且0<w<1,则SOR迭代法收敛。w=1时，SOR迭代法退化为Gauss-Seidel迭代法。

2 数值结果
2.1问题
考虑两点边值问题：

容易知道它的精确解为：
为了将微分方程离散，把[0,1]区间n等分，令h=1/n，，得到差分方程
，从而得到迭代方程组的系数矩阵A。
对=1，a=1/2,n=100,分别用jacobi，G-S，超松弛迭代法分别求线性方程组的解，要求4位有效数字，然后比较与精确解的误差。
对=0.1，=0.01，=0.001，考虑同样问题。

1.方程的表示及存储
由于本题中线性方程组的系数矩阵为三对角矩阵，所以可以采用紧缩方法存储，即

然后在矩阵乘法时对下标处理一下即可。但是考虑到三种迭代方法的一般性，且本题中n=200并不是很大，所以实验中并没有采用紧缩存储，而是采用了直接存储。
2. 边值条件的处理
由于差分得到的方程组的第一行和最后一行中分别出现了边值y(0)与y(1)作为常数项，因此要在常向量的第一项和最后一项作一些修改：

3.迭代终止条件
首先确定要求的精度tol ，我们希望当

则停止迭代。
对于迭代格式，若且，则迭代序列的
第k 次近似解和精确解之间有估计式。
由题目要求知我们需要有，而由上面的迭代估计，只要，即即可。而本题中q可近似取为，因此最后令迭代终止条件为

4.SOR 迭代中最佳松弛因子的选取
由于SOR 迭代法的效果和其松弛因子w的选取有关，所以有必要选取合适的松弛因子。当选择最佳松弛因子
时，SOR 方法的迭代速度最快。
Matlab实现：
迭代矩阵是n-1阶的，不是n阶；
等号右端向量b的最后一项，不是ah^2，而是ah^2-eps-h

2.2精确解

带入a=1/2,=1
代码：
>> clear
>> x=linspace(0,1);
truy=(1-0.5)/(1-exp(-1/1))*(1-exp(-x./1))+x.*0.5;
figure;
plot(x,truy,'g','LineWidth',1.5);
hold on;
Grid

图：

2.3三种迭代法
Jacobi法：代码见附录
Eps=1
结果：
迭代次数k：22273
结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

Eps=0.1
结果：
迭代次数k：8753
结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

Eps=0.01
结果：
迭代次数k：661
结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

G-S迭代法：代码见附录
Eps=1
结果：
迭代次数k：11125
结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

Eps=0.1
结果：
迭代次数k：4394
结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

Eps=0.01
结果：
迭代次数k：379
结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

超松弛法：代码见附录
Eps=1 w=1.56
结果：
迭代次数k：3503
结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

Eps=0.1 w=1.56
结果：
迭代次数k：1369
结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

Eps=0.01 w=1.56
结果：
迭代次数k：131
结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

3 分析讨论及心得体会
3.1三种方法的比较
Jacobi、G-S、超松弛法，三者都能够取得对精确解的良好逼近，但是，在相同的精度条件下，三者的收敛速度是不一样的，jacobi<G-S<超松弛，也就是说，在迭代次数相同的条件下，精度：jacobi<G-S<超松弛。
3.2心得体会
这次课程设计，平时感觉挺简单的那些枯燥单调的代码和数学公式，真正到了自己运用的时候却无从下手，但是，解决问题的过程恰是不断学习的过程：数学算法转换为代码的过程要对题目有深入的了解，然后对程序函数定义还要有一定的掌握能力，通过这个的过程让我巩固了自己的数学知识，对数学专业知识和MATLAB的操作有了更深的体会。 
 课程设计中遇到的问题只凭自己苦思冥想是不能全部解决的，这是同学老师的建议和网络给了我很大的帮助。遇到自己解决不了的问题时，多多向老师同学请教，或许问题就能迎刃而解。

4参考文献
[1]徐树方.数值线性代数.北京:北京大学出版社,1995.
[2]马昌凤.现代数值分析.北京:国防工业出版社.2013.
[3]刘春凤,米翠兰.实用数值分析教程.北京冶金工业出版社.2006

5附录
源代码
1.Jacobi：
function [y,k]=jacobi2(a,eps,h,delta)
n=1.0/h;
A=ones(n-1);
y=zeros(n-1,1);
z=zeros(n-1,1);
k=0;
for i=1:n-1
for j=1:n-1
A(i,j)=0;
end
end
for i=1:n-1
A(i,i)=-(2*eps+h);
end
for i=1:n-1
for j=1:n-1
if i==j+1
A(i,j)=eps;
end
if i==j-1
A(i,j)=eps+h;
end
end
end
b=zeros(n-1,1);
for i=1:n-2
b(i,1)=a*h^2;
end
b(n-1,1)=a*h^2-eps-h;
D=zeros(n-1);
for i=1:n-1
D(i,i)=A(i,i);
end
L=zeros(n-1);
for i=1:n-1
for j=1:n-1
if i>j
L(i,j)=-A(i,j);
end
end
end
U=zeros(n-1);
for i=1:n-1
for j=1:n-1
if i<j
U(i,j)=-A(i,j);
end
end
end
B=D\(L+U);
g=D\b;
while 1
z=B*y+g;
if norm(z-y,inf)<delta
break;
end
y=z;k=k+1;
end

x=linspace(0,1);
truy=(1-a)/(1-exp(-1/eps))*(1-exp(-x./eps))+x.*a;
figure;
plot(100*x,truy,'g','LineWidth',5);
hold on;
grid
hold on;
plot(y,'b')

2.G-S:
function [y,k]=gs2(a,eps,h,delta)
n=1.0/h;
A=ones(n-1);
y=zeros(n-1,1);
z=zeros(n-1,1);
k=0;
for i=1:n-1
for j=1:n-1
A(i,j)=0;
end
end
for i=1:n-1
A(i,i)=-(2*eps+h);
end
for i=1:n-1
for j=1:n-1
if i==j+1
A(i,j)=eps;
end
if i==j-1
A(i,j)=eps+h;
end
end
end
b=zeros(n-1,1);
for i=1:n-2
b(i,1)=a*h^2;
end
b(n-1,1)=a*h^2-eps-h;
D=zeros(n-1);
for i=1:n-1
D(i,i)=A(i,i);
end
L=zeros(n-1);
for i=1:n-1
for j=1:n-1
if i>j
L(i,j)=-A(i,j);
end
end
end
U=zeros(n-1);
for i=1:n-1
for j=1:n-1
if i<j
U(i,j)=-A(i,j);
end
end
end
B=D\(L+U);
g=D\b;
while 1
z=(D-L)\U*y+(D-L)\b;
if norm(z-y,inf)<delta
break;
end
y=z;k=k+1;
end

x=linspace(0,1);
truy=(1-a)/(1-exp(-1/eps))*(1-exp(-x./eps))+x.*a;
figure;
plot(100*x,truy,'g','LineWidth',5);
hold on;
grid
hold on;
plot(y,'b')

3.SOR：
function [y,k]=sor(a,eps,h,delta,w)
n=1.0/h;
A=ones(n-1);
y=zeros(n-1,1);
z=zeros(n-1,1);
k=0;
for i=1:n-1
for j=1:n-1
A(i,j)=0;
end
end
for i=1:n-1
A(i,i)=-(2*eps+h);
end
for i=1:n-1
for j=1:n-1
if i==j+1
A(i,j)=eps;
end
if i==j-1
A(i,j)=eps+h;
end
end
end
b=zeros(n-1,1);
for i=1:n-2
b(i,1)=a*h^2;
end
b(n-1,1)=a*h^2-eps-h;
D=zeros(n-1);
for i=1:n-1
D(i,i)=A(i,i);
end
L=zeros(n-1);
for i=1:n-1
for j=1:n-1
if i>j
L(i,j)=-A(i,j);
end
end
end
U=zeros(n-1);
for i=1:n-1
for j=1:n-1
if i<j
U(i,j)=-A(i,j);
end
end
end
B=D\(L+U);
g=D\b;
Lw=((D-w*L)^-1)*((1-w)*D+w*U);
while 1
z=Lw*y+w*(D-w*L)^-1*b;
if norm(z-y,inf)<delta
break;
end
y=z;k=k+1;
end

x=linspace(0,1);
truy=(1-a)/(1-exp(-1/eps))*(1-exp(-x./eps))+x.*a;
figure;
plot(100*x,truy,'g','LineWidth',5);
hold on;
grid
hold on;
plot(y,'b')