数值线性的代数课设课件资料文档格式.docx
《数值线性的代数课设课件资料文档格式.docx》由会员分享,可在线阅读,更多相关《数值线性的代数课设课件资料文档格式.docx(21页珍藏版)》请在冰豆网上搜索。
这三种方法的原理大致相同,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<
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);
holdon;
Grid
图:
2.3三种迭代法
Jacobi法:
代码见附录
Eps=1
结果:
迭代次数k:
22273
结果与精确解的比较图(绿色粗线是精确解,黑色细线是迭代结果)
Eps=0.1
8753
Eps=0.01
661
G-S迭代法:
11125
4394
379
超松弛法:
Eps=1w=1.56
3503
Eps=0.1w=1.56
1369
Eps=0.01w=1.56
131
3分析讨论及心得体会
3.1三种方法的比较
Jacobi、G-S、超松弛法,三者都能够取得对精确解的良好逼近,但是,在相同的精度条件下,三者的收敛速度是不一样的,jacobi<
G-S<
超松弛,也就是说,在迭代次数相同的条件下,精度:
jacobi<
超松弛。
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;
fori=1:
n-1
forj=1:
A(i,j)=0;
end
end
A(i,i)=-(2*eps+h);
ifi==j+1
A(i,j)=eps;
ifi==j-1
A(i,j)=eps+h;
b=zeros(n-1,1);
n-2
b(i,1)=a*h^2;
b(n-1,1)=a*h^2-eps-h;
D=zeros(n-1);
D(i,i)=A(i,i);
L=zeros(n-1);
ifi>
j
L(i,j)=-A(i,j);
U=zeros(n-1);
ifi<
U(i,j)=-A(i,j);
B=D\(L+U);
g=D\b;
while1
z=B*y+g;
ifnorm(z-y,inf)<
delta
break;
y=z;
k=k+1;
x=linspace(0,1);
truy=(1-a)/(1-exp(-1/eps))*(1-exp(-x./eps))+x.*a;
plot(100*x,truy,'
5);
grid
plot(y,'
b'
)
2.G-S:
function[y,k]=gs2(a,eps,h,delta)
z=(D-L)\U*y+(D-L)\b;
3.SOR:
function[y,k]=sor(a,eps,h,delta,w)