数值实验题三.docxVIP

3.用有限差分方法（五点差分格式）求解正方形域上的Poisson方程边值问题用MATLAB语言编写求解线性方程组的算法程序，采用下列方法计算，比较计算结果和算法性能，对计算结果给出结论。（1）用Jacobi迭代法求解线性方程组。（2）用块Jacobi迭代法求解线性方程组。（3）用（预条件）共轭斜向量法求解线性方程组。解：基本原理：由差分格式可得:写成矩阵形式：Au=f其中：其中：根据求得的结果分别采用普通Jacobi迭代法和块Jacobi迭代法分别求解，其中，块Jacobi迭代法还需要调用之前追赶法的迭代程序进行求解。具体程序如下：(1)普通Jacobi迭代法：function [u,k,t]=Jcb(n)%u 表示方程组的解； Ub 表示迭代矩阵；k 表示迭代次数；n 表示非边界点数% f 表示线性方程组A*U=f的右端矩阵f ；e 表示允许误差界；er 表示迭代误差% t 表示计算时间f(2:n+1,2:n+1)= 2*(n+1)^(-2); %初始化fUb=zeros(n+2); %初始化Ub为n+2阶零矩阵u=zeros(n+2); %初始化u为n+2阶零矩阵e=1e-9; %设置误差界tic;for k=1:1000 %迭代求解 er=0; for i=2:n+1 for j=2:n+1 Ub(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(Ub(i,j)-u(i,j)); %估计当前误差 end end u=Ub; %将矩阵Ub的值赋予u if er/n^2e,break; end %判断是否达到计算精度，如果达到则退出循环endt=toc; %获得计算时间输入命令：[u,k,t]=Jcb(9)，surf(u)结果：k=304 t=0.0026 u结果如下：（2）块Jacobi迭代法function [u,k,t]=KJcb(n)%块Jacobi迭代法% u 表示方程组的解；k 表示迭代次数；n 表示非边界点数；% f 表示线性方程组A*u=f的右端矩阵f； a 表示方程组系数矩阵的下对角线元素 % b 表示方程组系数矩阵的主对角线元素；c 表示方程组系数矩阵的上对角线元素；d 表示追赶法所求方程的右端向量% e 表示允许误差界；er 表示迭代误差；f(2:n+1,2:n+1)=2*(n+1)^(-2); %初始化fu=zeros(n+2); %初始化u为n+2阶零矩阵e=1e-9; %设置误差界 a=-1*ones(1,n);b=4*ones(1,n);c=-1*ones(1,n); tic;for k=1:2000 % 迭代求解 er=0; for j=2:n+1 Ub=u(:,j); d(1:n)=f(2:n+1,j)+u(2:n+1,j-1)+u(2:n+1,j+1); % 块Jacobi迭代 x=zg(a,b,c,d); %追赶法求解 er=er+norm(Ub-u(:,j),1); % 估计误差 end if er/n^2e ,break;end %判断是否达到计算精度，如果达到则退出循环endt=toc;追赶法程序采用之前编写的程序，如下：function x=zg(a,b,c,d)n=length(b);u(1)=b(1);for k=2:nif u(k-1)==0,D=0,return;end l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*c(k-1);end % 追赶法求解之追过程,求解Ly=dy(1)=d(1);for k=2:n y(k)=d(k)-l(k)*y(k-1);end % 追赶法求解之赶过程,求解Uz=yif u(n)==0 ,D=0,break; end x(n)=y(n) /u(n); endfor k=n-1:-1:1x(k)=(y(k)-c(k)*x(k+1))/u(k);end输入命令：[u,k,t]=KJcb(9)，surf(u)结果：k=86 t=0.0618 u结果如下： （3）共轭向量法function A=xishu(N) %存储离散化后非边界点构成的系数矩阵A=zeros(N^2);for i=1:N^2 A(i,i)=4; if i+NN^2+1 A(i,i+N)=-1; A(i+N,i)=A(i,i+N); end if mod(i,N)~=0 A(i,i+1)=-1; A(i+1,i)=A(i,i+1); endendfunctio