一、数值求解如下正方形域上的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,1i j i j i j i j i j ijj N j i i N u u u u u h f i j Nu 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|<ep 的绝对值，判断是否停机 6． 如果小于规定误差，迭代终止； 7． 输出结果u 和迭代次数k8． 在块的迭代中调用了追赶法的求解子程序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.011T T N N TN 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 IA -⎛⎫ ⎪-⎪= ⎪- ⎪-⎝⎭四、编写求解线性方程组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+1u(1,j)=(j-1)*h*(1-(j-1)*h);u(n+2,j)=(j-1)*h*(1-(j-1)*h);end%开始迭代while k<=it_maxv=u;for i=2:n+1for j=2:n+1u(i,j)=(v(i-1,j)+v(i+1,j)+v(i,j-1)+v(i,j+1)+f(i,j))/4;endendif max(abs(u-v))<epbreak;endk=k+1;end2．用块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=dendx(n)=y(n)/u(n); % 追赶法求解之赶过程求解Uz=yfor m=n-1:-1:1if u(m)==0 ,D=0,break; endx(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+1u(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+1d(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';endif max(abs(Ub-u))<epbreak;endk=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+1u(1,j)=(j-1)*h*(1-(j-1)*h);u(n+2,j)=(j-1)*h*(1-(j-1)*h);endwhile k<=it_maxuk=u;%用于存放的k次迭代的值for i=2:n+1for j=2:n+1u(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);endendif max(abs(uk-u))<epbreak;endk=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:nfor j=1:nif abs(i-j)==1AA(i,j)=1;endendendAB=eye(n);%安矩阵的块给A赋值for k=1:nfor i=1:nfor j=1:nA((i+(k-1)*n),(j+(k-1)*n))=AA(i,j);endendendfor k=1:nfor i=1:nfor j=1:nA((i+k*n),(j+(k-1)*n))=AB(i,j);A((i+(k-1)*n),(j+k*n))=AB(i,j);endendendfunction [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+1u(1,j)=(j-1)*h*(1-(j-1)*h);u(n+2,j)=(j-1)*h*(1-(j-1)*h);endwhile k<=it_maxuk=u;er=0;for i=1:nsum=zeros(n,1);for j=1:nsum=sum+A((((i-1)*n+1):i*n),(((j-1)*n+1):j*n))*u(2:n+1,j+1);endu(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);endif max(abs(er/n^2))<epbreak;endk=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+1u(1,j)=(j-1)*h*(1-(j-1)*h);u(n+2,j)=(j-1)*h*(1-(j-1)*h);endwhile k<=it_maxuk=u;%用于存放的k次迭代的值for i=2:n+1for j=2:n+1u(i,j)=(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)+f(i,j))/4;%用于存放的k+1次迭代的值endendif max(abs(u-uk))<epbreak;endk=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+1u(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+1d(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';endif max(abs(Ub-u))<epbreak;endk=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 00 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 00 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 00 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 00 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 00 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 00 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 00 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 00 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 00 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 00 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 k = 332tic;[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。