迭代法求两点边值问题课程设计分解.doc

迭代法求两点边值问题课程设计分解

28
2020年4月19日
文档仅供参考，不当之处，请联系改正。

迭代法求两点边值问题（一）

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

关键词：数值计算 迭代法 三种公式
ITERATION METHOD FOR TWO-POINT BOUNDARY VALUE PROBLEMS

ABSTRACT

Today's environment, numerical computation is more and more dependent on the computer. Many problems in the large-scale scientific computing and engineering solution, ultimately comes down to solving large sparse system of linear equations, the solution time occupies a large proportion in the whole problem solving time, some even reached 80%. Due to the current scientific research and large projects in various complex can claim to the computing precision and computing speed is higher and higher. Therefore, the linear algebraic equations, which is the basis for large-scale scientific computing efficient numerical solution has attracted widespread attention. To solve the general iterative method is used to this kind of equations. This time I do the iteration method of two-point boundary value problem is accomplished by MATLAB.
About the iteration method, there are many kinds of solution formula is: Jacobi, G - S and over relaxation iteration method. The principles of these three methods are roughly the same, Jacobi need at the beginning of a given vector, G - S requires the given initial value, a super relaxation method is weighted average of the transformation of Gauss’ - Seidel iterative method. And this article is for iteration solving large sparse linear equations with three kinds of iterative method (Jacobi, Gauss - Seidel and over relaxation iteration method) the error between the convergence speed and the exact solution of comparative research.

Key words: numerical calculation iterative method three kinds of formula

目 录
1 课程目的与要求……………………………………………………………………………1
1、1课程目标……………………………………………………………………………1
1、2 课程的实现环境……………………………………………………………………1
2 概要设计……………………………………………………………………………………2
3 详细设计……………………………………………………………………………………5
4 测试结果与心得体会………………………………………………………………………8
5 参考文献…………………………………………………………………………………15
附录…………………………………………………………………………………………16
源程序代码…………………………………………………………………………………16

1 课程的目的与要求

1、1课程目标
利用我们所学的数值分析的迭代法知识对线性方程组的问题进行求解，同时把迭代法的求解过程用MATLAB程序语言实现其算法功能。
对实验结果进行记录，对比与分析最后得出结论。

1、2课程的实现环境
硬件要求能运行Windows 9.X操作系统的微机系统、Matlab软件。

2概要设计
2、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迭代。

3详细设计
3.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
精确解：
带入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

图：

4测试结果与心得体会

4、1程序运行情况
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
结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

4、2三种方法的比较
Jacobi、G-S、超松弛法，三者都能够取得对精确解的良好逼近，可是，在相同的精度条件下，三者的收敛速度是不一样的，jacobi<G-S<超松弛，也就是说，在迭代次数相同的条件下，精度：jacobi<G-S<超松弛。
4、3心得体会
本次数值方法的题目是迭代法求两点边值问题。首先回去温习了数值方法的理论知识，把迭代法的原理弄清楚了，请教老师，上网查资料，忙的不亦乐乎。在设计的过程中我遇到的很多问题，在老师的帮助和自己的思考下还是很好的完成了。
事情就是如此，努力代表收获。这次课程设计，平时感觉挺简单的那些枯燥单调的代码和数学公式，真正到了自己运用的时候却无从下手，可是，解决问题的过程恰是不断学习的过程：数学算法转换为代码的过程要对题目有深入的了解，然后对程序函数定义还要有一定的掌握能力，经过这个的过程让我巩固了自己的数学知识，对数学专业知识和MATLAB的操作有了更深的体会。 
 课程设计中遇到的问题只凭自己苦思冥想是不能全部解决的，这是同学老师的建议和网络给了我很大的帮助。遇到自己解决不了的问题时，多多向老师同学请教，或许问题就能迎刃而解。另外，生在这个信息化社会，互联网的作用不言而喻，许多问题只要在网上搜索一下就能得到令人满意的解答。 

参考文献
[1] 王建卫．MATLAB 7.X 程序设计．北京：中国水利水电出版社， ．
[2] 李庆扬，王能超，易大义．数值分析[M] ．武汉：华中科技大学出版社， .7．
[3] 清华大学、北京大学计算方法编写组．计算方法[M] ．北京：科学出版社，1980．

附录
源程序代码：
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')

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')

超松弛：
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')