用SOR迭代法

一、数值求解如下正方形域上的Poisson 方程边值问

二、2222(,)2,0,1

(0,)(1,)(1),01(,0)(,1)0,

01u u f x y x y x y u y u y y y y u x u x x ???

??-+==<

==-≤≤??==≤≤?

二、用椭圆型第一边值问题的五点差分格式得到线性方程组为

2,1,1,,1,10,1,,0,141,?,?,?,?0,1

i j i j i j i j i j ij

j N j i i N u u u u u h f i j N

u u u u i j N -+-+++----=≤≤====≤≤+，

写成矩阵形式Au=f 。其中 三、基本原理

程序步骤：所有的步骤基本一致 1． 设置u ，n ，并给u ，n 赋初值； 2． While 语句循环，到的6步 3． Up 我第K 次迭代的值； 4． 分别进行计算，sum=0; 例如:

Jacobi :sun= sum+A(i,j)*Ub; SOR 和Gauss_Seidel= sum+A(i,j)*u; 各自进行相应的下不运算。

5． 计算|Up-u|

8． 在块的迭代中调用了追赶法的求解子程序zg ，在SOR 设计了A 得自动生成子程序

creat_matix.

1122N N v b v b u f v b ???? ? ? ?

?== ? ? ? ?????

4114114ii A -?? ?- ?= ?- ?

-?

? 11,12,1,121,22,2,21,2,,2211,12,1,121,22,2,221,2,,(,,...,),(,,...,),......,(,,...,)(,,...,)?,(,,...,)?,......,(,,...,)?

1,999,

0.10.01

1T T N N T

N N N N N T T N N T N N N N N v u u u v u u u v u u u b h f f f b h f f f b h f f f h N h N ====+=+=+===+取或则或,2,,1,2,...,i j f i j N

==

1122NN A I I A A I I

A -?? ?-

?= ?

- ?-?

?

四、编写求解线性方程组Au=f的算法程序，用下列方法编程计算，并比较计算速度。1．用Jacobi迭代法求解线性方程组Au=f。

function [u,k]=Jacobi(n,ep,it_max)

%JACOBI迭代法求Au=f;

%n迭代次数;

%ep为精度要求;

% it_max为最大迭代次数;

% u为方程组的解;

% k为迭代次数;

h=1/(n+1);

f(2:n+1,2:n+1)=h*h*2;%给f赋初值

u=zeros(n+2,n+2);v=zeros(n+2,n+2);k=1;

%给u赋初值

for j=2:n+1

u(1,j)=(j-1)*h*(1-(j-1)*h);

u(n+2,j)=(j-1)*h*(1-(j-1)*h);

end

%开始迭代

while k<=it_max

v=u;

for i=2:n+1

for j=2:n+1

u(i,j)=(v(i-1,j)+v(i+1,j)+v(i,j-1)+v(i,j+1)+f(i,j))/4;

end

end

if max(abs(u-v))

break;

end

k=k+1;

end

2．用块Jacobi迭代法求解线性方程组Au=f。

function x=zg(a,b,c,d)

%求解三对角方程的追赶法

n=length(b);

u(1)=b(1);y(1)=d(1);

for i=2:n l(i)=a(i)/u(i-1);

u(i)=b(i)-l(i)*c(i-1);

y(i)=d(i)-l(i)*y(i-1); % 追赶法求解之追过程求解Ly=d

end

x(n)=y(n)/u(n); % 追赶法求解之赶过程求解Uz=y

for m=n-1:-1:1

if u(m)==0 ,D=0,break; end

x(m)=(y(m)-c(m)*x(m+1))/u(m);

end

function [u,k]=Jacobi_block(n,ep,it_max)

% 用块jacobi迭代法求解线性方程组A*u=f

% u: 方程组的解；k: 迭代次数；n: 非边界点数；

% a: 方程组系数矩阵的下对角线元素；b: 方程组系数矩阵的主对角线元素；

% c: 方程组系数矩阵的上对角线元素；d: 追赶法所求方程的右端向量；% ep: 允许误差界; %it_max:迭代的最大次数；% function x=zg(a,b,c,d) 子函数追赶法求解；

h=1/(n+1);

f(2:n+1,2:n+1)=h*h*2;

a=-1*ones(1,n); b=4*ones(1,n);c=-1*ones(1,n);k=1;

u=zeros(n+2,n+2);

%给u赋初值

for j=2:n+1

u(1,j)=(j-1)*h*(1-(j-1)*h);

u(n+2,j)=(j-1)*h*(1-(j-1)*h);

end

%开始迭代

while k<=it_max

Ub=u;

for j=2:n+1

d(1:n)=f(2:n+1,j)+Ub(2:n+1,j-1)+Ub(2:n+1,j+1) ;

x=zg(a,b,c,d); % 调用子函数追赶法求解

u(2:n+1,j)=x';

end

if max(abs(Ub-u))

break;

end

k=k+1;

end

3．用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。

function [u,k]=SOR(n,ep,w,it_max)

%SOR迭代法%n迭代次数%ep为精度要求% it_max为最大迭代次数

% u为方程组的解% k为迭代次数%w为松弛因子

h=1/(n+1);

f(2:n+1,2:n+1)=h*h*2;

u=zeros(n+2,n+2);k=1;

for j=2:n+1

u(1,j)=(j-1)*h*(1-(j-1)*h);

u(n+2,j)=(j-1)*h*(1-(j-1)*h);

end

while k<=it_max

uk=u;%用于存放的k次迭代的值

for i=2:n+1

for j=2:n+1

u(i,j)=w*(-4*u(i,j)+u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)+f(i,j))/4;

u(i,j)=u(i,j)+uk(i,j);

end

end

if max(abs(uk-u))

break;

end

k=k+1;

end

4．用块SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。function [AA,A]=creat_matrix(n)

%自动的生成矩阵A

A=zeros((n)^2,(n)^2);%定义A的对角的对成NXN的方阵

AA=4*eye(n);

for i=1:n

for j=1:n

if abs(i-j)==1

AA(i,j)=1;

end

end

end

AB=eye(n);

%安矩阵的块给A赋值

for k=1:n

for i=1:n

for j=1:n

A((i+(k-1)*n),(j+(k-1)*n))=AA(i,j);

end

end

end

for k=1:n

for i=1:n

for j=1:n

A((i+k*n),(j+(k-1)*n))=AB(i,j);

A((i+(k-1)*n),(j+k*n))=AB(i,j);

end

end

end

function [u,k]=SOR_block(n,w,ep,it_max)

% 用块SOR迭代法求解线性方程组A*u=f

% u: 方程组的解；k: 迭代次数；n: 非边界点数；% ep: 允许误差界;

%it_max:迭代的最大次数；%A=creat-matrix（n），为创建A矩阵的子函数；[AA,A]=creat_matrix(n); %调用子函数；

h=1/(n+1);

k=1;

f(2:n+1,2:n+1)=h*h*2;

u=zeros(n+2,n+2);

for j=2:n+1

u(1,j)=(j-1)*h*(1-(j-1)*h);

u(n+2,j)=(j-1)*h*(1-(j-1)*h);

end

while k<=it_max

uk=u;

er=0;

for i=1:n

sum=zeros(n,1);

for j=1:n

sum=sum+A((((i-1)*n+1):i*n),(((j-1)*n+1):j*n))*u(2:n+1,j+1);

end

u(2:n+1,i+1)=uk(2:n+1,i+1)+w*inv(AA)*(f(2:n+1,i+1)-sum);

er=er+norm(uk(:,i+1)-u(:,i+1),1);

end

if max(abs(er/n^2))

break;

end

k=k+1;

end

5．用Gauss-Seidel迭代法求解线性方程组Au=f。

function [u,k]=Gauss_Seidel(n,ep,it_max)

%G--s迭代法%n迭代次数%ep为精度要求% it_max为最大迭代次数

% u为方程组的解% k为迭代次数

h=1/(n+1);

f(2:n+1,2:n+1)=h*h*2;

u=zeros(n+2,n+2);k=1;

for j=2:n+1

u(1,j)=(j-1)*h*(1-(j-1)*h);

u(n+2,j)=(j-1)*h*(1-(j-1)*h);

end

while k<=it_max

uk=u;%用于存放的k次迭代的值

for i=2:n+1

for j=2:n+1

u(i,j)=(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)+f(i,j))/4;%用于存放的k+1次迭代的值end

end

if max(abs(u-uk))

break;

end

k=k+1;

end

6．用块Gauss-Seidel迭代法求解线性方程组Au=f。

function [u,k]=Gauss_Seidel_block(n,ep,it_max)

% 用块Gauss-seidel迭代法求解线性方程组A*u=f

% u: 方程组的解；k: 迭代次数；n: 非边界点数；

% a: 方程组系数矩阵的下对角线元素；b: 方程组系数矩阵的主对角线元素；

% c: 方程组系数矩阵的上对角线元素；d: 追赶法所求方程的右端向量；% ep: 允许误差界; %it_max:迭代的最大次数；% function x=zg(a,b,c,d) 子函数追赶法求解；

h=1/(n+1);

f(2:n+1,2:n+1)=h*h*2;

a=-1*ones(1,n); b=4*ones(1,n);c=-1*ones(1,n);k=1;

u=zeros(n+2,n+2);

for j=2:n+1

u(1,j)=(j-1)*h*(1-(j-1)*h);

u(n+2,j)=(j-1)*h*(1-(j-1)*h);

end

% for j=2:n+1

% f(2,j)=h*h+(j-1)*h*(1-(j-1)*h);

% f(n+1,j)=h*h+(j-1)*h*(1-(j-1)*h);

% end

while k<=it_max

Ub=u;

for j=2:n+1

d(1:n)=f(2:n+1,j)+u(2:n+1,j-1)+u(2:n+1,j+1) ;

x=zg(a,b,c,d); % 用追赶法求解

u(2:n+1,j)=x';

end

if max(abs(Ub-u))

break;

end

k=k+1;

end

五、各种算法的实验结果对比

在MA TLAB中输入

n=9; ep=0.000000001; it_max=1000; w=1.8;

1)[u,k]=Jacobi(n,ep,it_max) 回车

u =

0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0

0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0

0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0

0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0

0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0

0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0

0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0

0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0

0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0

0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0

0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 k = 332

tic;[u,k]= (n,ep,it_max);toc;

Elapsed time is 0.011470 seconds.

2)[u,k]=Gauss_Seidel(n,ep,it_max) 回车

k = 174

>> tic;[u,k]= (n,ep,it_max);toc;

Elapsed time is 0.006760 seconds.

3)[u,k]= (n,ep,w,it_max) 回车

k =91

>> tic;[u,k]=SOR(n,w,ep,it_max);toc;

Elapsed time is 0.000497 seconds.

把松弛系数w调整为0.8,1.0,1.1, 1,7,1.8发现；迭代次数在w=1。8时最少，为91步。

4) [u,k]=Jacobi_block(n,ep,it_max)

u =

0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0

0 0.0256 0.0413 0.0508 0.0560 0.0577 0.0560 0.0508 0.0413 0.0256 0

0 0.0413 0.0686 0.0859 0.0955 0.0986 0.0955 0.0859 0.0686 0.0413 0

0 0.0508 0.0859 0.1088 0.1216 0.1258 0.1216 0.1088 0.0859 0.0508 0

0 0.0560 0.0955 0.1216 0.1364 0.1412 0.1364 0.1216 0.0955 0.0560 0

0 0.0577 0.0986 0.1258 0.1412 0.1462 0.1412 0.1258 0.0986 0.0577 0

0 0.0560 0.0955 0.1216 0.1364 0.1412 0.1364 0.1216 0.0955 0.0560 0

0 0.0508 0.0859 0.1088 0.1216 0.1258 0.1216 0.1088 0.0859 0.0508 0

0 0.0413 0.0686 0.0859 0.0955 0.0986 0.0955 0.0859 0.0686 0.0413 0

0 0.0256 0.0413 0.0508 0.0560 0.0577 0.0560 0.0508 0.0413 0.0256 0

0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0

k = 170

>> tic;[u,k]=Jacobi_block(n,ep,it_max);toc;

Elapsed time is 0.159638 seconds.

5)[u,k]=Gauss_Seidel_block(n,ep,it_max)

k =90

>> tic;[u,k]=Gauss_Seidel_block(n,ep,it_max);toc;

Elapsed time is 0.078421 seconds.

6)[u,k]=SOR_block(n,w,ep,it_max)

k =71

>> tic;[u,k]=SOR_block(n,w,ep,it_max);toc;

Elapsed time is 0.098572 seconds.

把松弛系数w调整为0.8,1.0,1.1,1,2 1,3发现；迭代次数在w=11时最少，为33步。

分析结果：

1．块的迭代的生成最终结果相同，一般迭代的最终结果相同，这可能是由于误差判断的方式决定的，快的迭代为均差，而一般迭代为单个元素最大的差值。

2．由上面1-3的结果课已看出，Jacobi, Gauss_Seidel,SOR在相同的条件下，迭代的次数依次递减；耗时依次递减。

3．由上面4-6的结果课已看出，用块的Jacobi, Gauss_Seidel,SOR在相同的条件下，迭代的次数依次递减；块SOR的耗时最长。应为SOR块的迭代是构造AX=b的方

式实现的，没法和前两中块的迭代作比较。

4．W对SOR迭代还是对块的SOR迭代的影响很大，只有找到各自合适的w时才能达到最少迭代次数和耗时。