第三章编程作业

1、数值实验题三3.用有限差分方法（五点差分格式）求解正方形域上的Poisson方程边值问题用MATLAB语言编写求解线性方程组的算法程序，采用下列方法计算，比较计算结果和算法性能，对计算结果给出结论。（1）用Jacobi迭代法求解线性方程组。（2）用块Jacobi迭代法求解线性方程组。（3）用（预条件）共轭斜向量法求解线性方程组。解：由差分格式可得:写成矩阵形式：Au=f其中：其中：（1）Jacobi迭代法：function U,k,er,t=Jcb(n)% Jacobi迭代法% U 表示方程组的解；h 表示步长； A 表示迭代矩阵；k 表示迭代次数；n 表示非边界点数% f 表示线性方程组A*

2、U=f的右端矩阵f ；e 表示允许误差界；er 表示迭代误差% t 表示计算时间h=1/(n+1);f(2:n+1,2:n+1)= 2*h2; %初始化fA=zeros(n+2); %初始化A为n+2阶零矩阵U=zeros(n+2); %初始化U为n+2阶零矩阵e=0.; %设置误差界tic;for k=1:10000 %迭代求解 er=0; for i=2:n+1 for j=2:n+1 A(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); %估计当前误差 end end U=A; %将矩

3、阵A的值赋予U er=er/n2; if er 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

4、 表示系数矩阵A所分解成的下三角阵L中的下对角线元素 l(i)；z 表示系数矩阵A所分解成的上三角阵U中的主对角线元素 z(i)h=1/(n+1);f(2:n+1,2:n+1)=2*h2; %初始化fA=zeros(n+2); %初始化A为n+2阶零矩阵U=zeros(n+2); %初始化U为n+2阶零矩阵e=0.; %设置误差界 a=-1*ones(1,n);b=4*ones(1,n);c=-1*ones(1,n); tic;for k=1:1000 % 迭代求解 er=0; for j=2:n+1 q=zeros(1,n); d=f(2:n+1,j)+U(2:n+1,j-1)+U(2:n+

5、1,j+1); % 块Jacobi迭代 d=d; A(2:n+1,j)=zg(a,b,c,d); er=er+norm(A(:,j)-U(:,j),1); % 估计误差 end U=A; er=er/n2; if ere ,break;end %判断是否达到计算精度，如果达到则退出循环endt=toc;追赶法：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

6、=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运行结果：U,k,er,t=KJcb(9)其中U如下：k = 163er = 9.5903e-010t = 0.1329（3）共轭斜量法function A=xishu(N) %存储离散化后非边界点构成的系数矩阵A=zeros(N2);for i=1:N2 A(i,i)=4; if i+N1e-10 w=a*p; t=q0/(p*w); x=x+t*p; r=r-t*

7、w; q=r*r; s=q/q0; p=r+s*p; q0=q;end运行： a=xishu(9);for i=1:92b(i,1)=0.02;x(i,1)=0;end x=cg3(a,b,x)reshape(x,9,9)结果：下面的数据均为非边界点的值。4.用有限差分方法（五点差分格式）求解正方形域上的Poisson方程边界问题用Matlab语言编写求解线性方程组的算法程序，采用下列方法计算，并比较方法的计算速度。（1）用SOR迭代法求解线性方程组，用试算法确定最佳松弛因子。（2）用块SOR迭代法求解线性方程组，用试算法确定最佳松弛因子。（3）用（预条件）共轭斜向量法求解差分方程组。解：（1

8、）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*h2; %初始化fU=zeros(n+2); %初始化U矩阵for l=1:n+2 U(1,l)=(l-1)*h*(1-(l-1)*h); U(n+2,l)=(l-1)*h*(1-(l-1)*h);ende=0.; %设置误差界 tic;for k=1:10000 %迭代求解 er=0

9、; for i=2:n+1 for j=2:n+1 Ub=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); %估计当前误差 end end er=er/n2; if ere,break; %判断是否达到计算精度，如果达到则退出循环 endendt=toc; %获得计算时间最佳松弛因子：function goodw,goodk=zuijiaw(n)%计算最佳w%精度为0.01U,k,er,t=SOR(n,0.1);goodk=k;goodw=0.01

10、;for w=0.02:0.01:2 U,k,er,t=SOR(n,w); if k=goodk goodk=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所分解成的下三

11、角阵L中的下对角线元素 l(i)；z% 表示系数矩阵A所分解成的上三角阵U中的主对角线元素 z(i)h=1/(n+1);f(2:n+1,2:n+1)=3*h2; %初始化fU=zeros(n+2); %初始化U矩阵for l=1:n+2 U(1,l)=(l-1)*h*(1-(l-1)*h); U(n+2,l)=(l-1)*h*(1-(l-1)*h);ende=0.; %设置误差界 a=-1*ones(1,n);b=4*ones(1,n);c=-1*ones(1,n); tic;for k=1:1000 % 迭代求解 er=0; for j=2:n+1 Ub=U(:,j); q=zeros(1,

12、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); for i=1:n U(i+1,j)=w*x(i)+(1-w)*U(i+1,j); end er=er+norm(Ub-U(:,j),1); % 估计误差 end er=er/n2; if ere ,break;end %判断是否达到计算精度，如果达到则退出循环endt=toc;最佳松弛因子：function goodw,goodk=zuijiaw(n)%计算最佳w%精度为0.01U,k,er,t=KSOR(n,0.01);goodk=k;goodw=0.01;for w=0.02:0.01:2 U,k,er,t=KSOR(n,w); if k=goodk goodk=k; goodw=w; endend（3）共轭斜量法