用SOR迭代法

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

2u2u22f(x,y)2,0x,y1yx二、 u(0,y)u(1,y)y(1y),0y10x1u(x,0)u(x,1)0,

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

4ui,jui1,jui1,jui,j1ui,j1h2fij

u0,j?,uN1,j?,ui,01i,jN?,ui,N1?，

0i,jN1写成矩阵形式Au=f。其中 A11IAIA22I

IANNv1vu2vNb1bf2bN4114Aii114v1(u1,1,u2,1,...,uN,1)T,v2(u1,2,u2,2,...,uN,2)T,......,vN(u1,N,u2,N,...,uN,N)Tb1h2(f1,1,f2,1,...,fN,1)T?,b2h2(f1,2,f2,2,...,fN,2)T?,......,bNh2(f1,N,f2,N,...,fN,N)T?1h,取N9或99,N1fi,j2,i,j1,2,...,N三、基本原理

程序步骤：所有的步骤基本一致 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.

则h0.1或0.01四 、编写求解线性方程组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))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))用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))用块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))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))用块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))五、各种算法的实验结果对比 在MATLAB中输入

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 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600

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;

0.0900 0 0.0900 0 0.0900 0 0.0900 0 0.0900 0 0.0900 0 0.0900 0 0.0900 0 0.0900 0 0.0900 0 0.0900 0 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时才能达到

最少迭代次数和耗时。