第三章编程作业

数值实验题三3.用有限差分方法（五点差分格式）求解正方形域上的Poisson方程边值问题2222(,)2,0,1(0,)(1,)(,0)(,1)0uufxyxyxyuyuyuxux用MATLAB语言编写求解线性方程组Auf的算法程序，采用下列方法计算，比较计算结果和算法性能，对计算结果给出结论。（1）用Jacobi迭代法求解线性方程组Auf。（2）用块Jacobi迭代法求解线性方程组Auf。（3）用（预条件）共轭斜向量法求解线性方程组Auf。解：由差分格式可得:2,1,1,,1,10,1,,0,14(1,)0(0,)ijijijijijijjNjiiNuuuuuhfijNuuuuijN，写成矩阵形式：Au=f11112222......,,...........NNNNAIvbIAvbAvbIIAvb其中：4114....,.....114NNiiiiAAR其中：11,12,1,121,22,2,21,2,,2211,12,1,121,22,2,221,2,,(,,...,),(,,...,),......,(,,...,)(,,...,),(,,...,),......,(,,...,)TTNNTNNNNNTTNNTNNNNNvuuuvuuuvuuubhfffbhfffbhfff（1）Jacobi迭代法：function[U,k,er,t]=Jcb(n)%Jacobi迭代法%U表示方程组的解；h表示步长；A表示迭代矩阵；k表示迭代次数；n表示非边界点数%f表示线性方程组A*U=f的右端矩阵f；e表示允许误差界；er表示迭代误差%t表示计算时间h=1/(n+1);f(2:n+1,2:n+1)=2*h^2;%初始化fA=zeros(n+2);%初始化A为n+2阶零矩阵U=zeros(n+2);%初始化U为n+2阶零矩阵e=0.000000001;%设置误差界tic;fork=1:10000%迭代求解er=0;fori=2:n+1forj=2:n+1A(i,j)=(U(i-1,j)+U(i+1,j)+U(i,j+1)+U(i,j-1)+f(i,j))/4;er=er+abs(A(i,j)-U(i,j));%估计当前误差endendU=A;%将矩阵A的值赋予Uer=er/n^2;ifere,break;%判断是否达到计算精度，如果达到则退出循环endendt=toc;%获得计算时间[U,k,er,t]=Jcb(9)k=304er=9.7771e-010t=0.0077U如下：对U作图：（2）块Jacobi方法：function[U,k,er,t]=KJcb(n)%块Jacobi迭代法%u表示方程组的解h表示步长；k表示迭代次数；n表示非边界点数；%f表示线性方程组A*u=f的右端矩阵f；q表示n+2维向量；a表示方程组系数矩阵的下对角线元素%b表示方程组系数矩阵的主对角线元素；c表示方程组系数矩阵的上对角线元素；d表示追赶法所求方程的右端向量%e表示允许误差界；er表示迭代误差；l表示系数矩阵A所分解成的下三角阵L中的下对角线元素l(i)；z表示系数矩阵A所分解成的上三角阵U中的主对角线元素z(i)h=1/(n+1);f(2:n+1,2:n+1)=2*h^2;%初始化fA=zeros(n+2);%初始化A为n+2阶零矩阵U=zeros(n+2);%初始化U为n+2阶零矩阵e=0.000000001;%设置误差界a=-1*ones(1,n);b=4*ones(1,n);c=-1*ones(1,n);tic;fork=1:1000%迭代求解er=0;forj=2:n+1q=zeros(1,n);d=f(2:n+1,j)+U(2:n+1,j-1)+U(2:n+1,j+1);%块Jacobi迭代d=d';A(2:n+1,j)=zg(a,b,c,d);er=er+norm(A(:,j)-U(:,j),1);%估计误差endU=A;er=er/n^2;ifere,break;end%判断是否达到计算精度，如果达到则退出循环endt=toc;追赶法：functionx=zg(a,b,c,d)n=length(b);u(1)=b(1);y(1)=d(1);fori=2:nl(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=yform=n-1:-1:1ifu(m)==0,D=0,break;endx(m)=(y(m)-c(m)*x(m+1))/u(m);end运行结果：[U,k,er,t]=KJcb(9)其中U如下：k=163er=9.5903e-010t=0.1329（3）共轭斜量法functionA=xishu(N)%存储离散化后非边界点构成的系数矩阵A=zeros(N^2);fori=1:N^2A(i,i)=4;ifi+NN^2+1A(i,i+N)=-1;A(i+N,i)=A(i,i+N);endifmod(i,N)~=0A(i,i+1)=-1;A(i+1,i)=A(i,i+1);endendfunctionx=cg(a,b,x)%共轭向量法求解线性方程组%初始输入变量a：系数矩阵b：方程右端向量x输入的初始解n=length(b);r=b-a*x;p=r;q0=r'*r;whileq01e-10w=a*p;t=q0/(p'*w);x=x+t*p;r=r-t*w;q=r'*r;s=q/q0;p=r+s*p;q0=q;end运行：a=xishu(9);fori=1:9^2b(i,1)=0.02;x(i,1)=0;endx=cg3(a,b,x)reshape(x,9,9)结果：下面的数据均为非边界点的值。4.用有限差分方法（五点差分格式）求解正方形域上的Poisson方程边界问题2222()(,)3(0,1)(0,)(1,)0,(,0)(,1)(1)uufxyxyxyuyuyuxuxxx用Matlab语言编写求解线性方程组Aub的算法程序，采用下列方法计算，并比较方法的计算速度。（1）用SOR迭代法求解线性方程组Aub，用试算法确定最佳松弛因子。（2）用块SOR迭代法求解线性方程组Aub，用试算法确定最佳松弛因子。（3）用（预条件）共轭斜向量法求解差分方程组Aub。解：（1）SORfunction[U,k,er,t]=SOR(n,w)%SOR迭代法%U表示方程组的解；h表示步长；k表示迭代次数；n表示非边界点数%f表示线性方程组A*U=f的右端矩阵f；e表示允许误差界；er表示迭代误差%t表示计算时间h=1/(n+1);f(2:n+1,2:n+1)=3*h^2;%初始化fU=zeros(n+2);%初始化U矩阵forl=1:n+2U(1,l)=(l-1)*h*(1-(l-1)*h);U(n+2,l)=(l-1)*h*(1-(l-1)*h);ende=0.000000001;%设置误差界tic;fork=1:10000%迭代求解er=0;fori=2:n+1forj=2:n+1Ub=U(i,j);U(i,j)=w*((U(i-1,j)+U(i+1,j)+U(i,j+1)+U(i,j-1)+f(i,j))/4)+(1-w)*U(i,j);er=er+abs(Ub-U(i,j));%估计当前误差endender=er/n^2;ifere,break;%判断是否达到计算精度，如果达到则退出循环endendt=toc;%获得计算时间最佳松弛因子：function[goodw,goodk]=zuijiaw(n)%计算最佳w%精度为0.01[U,k,er,t]=SOR(n,0.1);goodk=k;goodw=0.01;forw=0.02:0.01:2[U,k,er,t]=SOR(n,w);ifk=goodkgoodk=k;goodw=w;endend（2）块SORfunction[U,k,er,t]=KSOR(n,w)%块SOR迭代法%u表示方程组的解h表示步长；k表示迭代次数；n表示非边界点数；%f表示线性方程组A*u=f的右端矩阵f；q表示n+2维向量；a表示方程组系数矩阵的下对角线元素%b表示方程组系数矩阵的主对角线元素；c表示方程组系数矩阵的上对角线元素；d表示追赶法所求方程的右端向量%e表示允许误差界；er表示迭代误差；l表示系数矩阵A所分解成的下三角阵L中的下对角线元素l(i)；z%表示系数矩阵A所分解成的上三角阵U中的主对角线元素z(i)h=1/(n+1);f(2:n+1,2:n+1)=3*h^2;%初始化fU=zeros(n+2);%初始化U矩阵forl=1:n+2U(1,l)=(l-1)*h*(1-(l-1)*h);U(n+2,l)=(l-1)*h*(1-(l-1)*h);ende=0.000000001;%设置误差界a=-1*ones(1,n);b=4*ones(1,n);c=-1*ones(1,n);tic;fork=1:1000%迭代求解er=0;forj=2:n+1Ub=U(:,j);q=zeros(1,n);d=f(2:n+1,j)+U(2:n+1,j-1)+U(2:n+1,j+1);%块SOR迭代d=d';x=zg(a,b,c,d);fori=1:nU(i+1,j)=w*x(i)+(1-w)*U(i+1,j);ender=er+norm(Ub-U(:,j),1);%估计误差ender=er/n^2;ifere,break;end%判断是否达到计算精度，如果达到则退出循环endt=toc;最佳松弛因子：function[goodw,goodk]=zuijiaw(n)%计算最佳w%精度为0.01[U,k,er,t]=KSOR(n,0.01);goodk=k;goodw=0.01;forw=0.02:0.01:2[U,k,er,t]=KSOR(n,w);ifk=goodkgoodk=k;goodw=w;endend（3）共轭斜量法functionA=xishu(N)%存储离散化后非边界点构成的系数矩阵A=zeros(N^2);fori=1:N^2A(i,i)=4;ifi+NN^2+1A(i,i+N)=-1;A(i+N,i)=A(i,i+N);endifmod(i,N)~=0A(i,i+1)=-1;A(i+1,i)=A(i,i+1);endendfunctionU=cg(n)%共轭向量法求解线性方程组h=1/(n+1);U=zeros(n+2);%初始化U矩阵forl=1:n+2U(1,l)=(l-1)*h*(1-(l-1)*h);U(n+2,l)=(l-1)*h*(1-(l-1)*h);endfori=1:n^2b(i,1)=0.03;x(i,1)=0;enda=xishu(n);n=length(b);r=b-a*x;p=r;q0=r'*r;whileq01e-10w=a*p;t=q0/(p'*w);x=x+t*p;r=r-t*w;q=r'*r;s=q/q0;p=r+s*p;q0=q;endN=n^0.5xx=reshape(x,N,N);U(2:N+1,2:N+1)=xx;