用SOR迭代法.docVIP

数值求解如下正方形域上的Poisson方程边值问 二、用椭圆型第一边值问题的五点差分格式得到线性方程组为 写成矩阵形式Au=f。其中 三、基本原理 程序步骤：所有的步骤基本一致 设置u，n，并给u，n赋初值； While 语句循环，到的6步 Up我第K次迭代的值； 分别进行计算，sum=0; 例如: Jacobi :sun= sum+A(i,j)*Ub; SOR和Gauss_Seidel= sum+A(i,j)*u; 各自进行相应的下不运算。 计算|Up-u|ep的绝对值，判断是否停机 如果小于规定误差，迭代终止； 输出结果u和迭代次数k 在块的迭代中调用了追赶法的求解子程序zg，在SOR设计了A得自动生成子程序creat_matix. 四 、编写求解线性方程组Au=f的算法程序，用下列方法编程计算，并比较计算速度。 用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 用块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))ep break; end k=k+1; end 用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。 function [u,k]=SOR(n,ep,w,it_max) %SOR迭代法 %n迭代次数 %ep为精度