用SOR迭代法.doc

一、 数值求解如下正方形域上的Poisson方程边值问 二、二、用椭圆型第一边值问题的五点差分格式得到线性方程组为写成矩阵形式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|ep的绝对值，判断是否停机6 如果小于规定误差，迭代终止；7 输出结果u和迭代次数k8 在块的迭代中调用了追赶法的求解子程序zg，在SOR设计了A得自动生成子程序creat_matix.四 、编写求解线性方程组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)ep 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); endfunction 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_maxUb=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)ep break; end k=k+1;end3 用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);endwhile 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)ep break; end k=k+1;end4 用块SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。function AA,A=creat_matrix(n)%自动的生成矩阵AA=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 endendAB=eye(n);%安矩阵的块给A赋值for k=1:nfor i=1:n for j=1:n A(i+(k-1)*n),(j+(k-1)*n)=AA(i,j); endendendfor 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 endendfunction 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);endwhile 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/n2)ep break; end k=k+1;end5 用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);endwhile 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)ep break; end k=k+1;end6 用块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);% endwhile k=it_maxUb=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) 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 0k = 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 块的