数值上机作业.docx

数值实验题一班级：研开发14-1班 姓名：马高强 学号：20142120904.用有限差分方法（五点差分格式）求解正方形域上的Poisson方程边值问题 用MATLAB语言编写求解线性方程组的算法程序，采用下列方法计算，比较计算结果和算法性能，对计算结果给出结论。（1）用SOR迭代法求解线性方程组，用试算法确定最佳松弛因子。（2）用块SOR迭代法求解线性方程组，用试算法确定最佳松弛因子。（3）用（预条件）共轭斜向量法求解线性方程组。解：首先，剖分求解区域。用平行于坐标轴的直线然后用数值微分公式对微分算子进行离散。既有其中表示的近似值。得到在每个点上的有限差分方程为：又称为五点差分格式.其中。在边界上有可得差分格式为：将差分方程组写成矩阵形式为，其中。如果把每一条线上的网格点看作一个组，如对矩阵也相应地分块其中则差分方程等价于。这样就可以用块迭代法进行求解。（1）用SOR迭代法求解线性方程组，用试算法确定最佳松弛因子。SOR迭代法：function u,k,t=SOR(n,w)%用SOR迭代法求解Poisson方程%输出结果：u，k，t 输入：n，wh=1/(n+1);f(2:n+1,2:n+1)=3*h2;%定义右端项u=zeros(n+2,n+2); %对边界进行差分 for i=1:n+2 u(i,1)=(i-1)*h*(1-(i-1)*h); u(i,n+2)=(i-1)*h*(1-(i-1)*h); ende=0.000000001; tic%程序开始的时间 %SOR迭代法过程for k=1:10000 er=0; 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 if er/n2e break; endendt1=toc%程序结束的时间t=max(t1);试算法确定最佳松弛因子：function goodw,goodk=Bestw(n)%试算法确定最佳松弛因子 u,k=SOR(n,0.01);goodw=0.01; goodk=k; for w=0.01:0.001:2 u,k=SOR(n,w); if k goodw,goodk=Bestw(9)goodw = 1.5330goodk =33 u,k,t=SOR(9,goodw)k = 33t = 2.3291e-004u的结果如下表：将u画成三维图 surf(u)（2）用块SOR迭代法求解线性方程组，用试算法确定最佳松弛因子。块SOR迭代法：function u,k,t=BSOR(n,w) %用块SOR求解Poisson问题，其中用到追赶法%输出结果：u，k，t 输入：n，wh=1/(n+1);f(2:n+1,2:n+1)=3*h2;u=zeros(n+2,n+2);a=-1*ones(1,n);b=4*ones(1,n);c=-1*ones(1,n); for i=1:n+2%对边界条件进行差分 u(i,1)=(i-1)*h*(1-(i-1)*h); u(i,n+2)=(i-1)*h*(1-(i-1)*h); end tic;%块SOR迭代法过程e=0.000000001; for k=1:1000 er=0; for j=2:n+1 ub=u(:,j); d=f(2:n+1,j)+u(2:n+1,j-1)+u(2:n+1,j+1); 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 if er/n2e break; endendt1=toc;t=max(t1);追赶法：function x=zg(a,b,c,d)n=length(b);%LU分解u(1)=b(1);for k=2:n if 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%求解Ux=yif u(n)=0 D=0; return;endx(n)=y(n)/u(n);for k=n-1:-1:1 x(k)=(y(k)-c(k)*x(k+1)/u(k);end试算法确定最佳松弛因子：function goodw,goodk=BBestw(n)%确定块SOR的最佳松弛因子 u,k=BSOR(n,0.01);goodw=0.01; goodk=k; for w=0.01:0.01:2 u,k=BSOR(n,w); if k goodw,goodk=BBestw(9)goodw = 1.4200goodk = 25 u,k,t=BSOR(9,goodw)k = 25t =0.0156u的结果如下表：将u画成三维图 surf(u)（3）用（预条件）共轭斜向量法求解线性方程组。共轭斜量法：function u,k,t=CG(n) %共轭斜量法求解Poisson问题A=zeros(n2,n2);for i=1:n2%系数 A(i,i)=4; if i+nn2+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); endend h=1/(n+1);u=zeros(n+2,n+2); for i=1:n+2%边界条件的差分 u(i,1)=(i-1)*h*(1-(i-1)*h); u(i,n+2)=(i-1)*h*(1-(i-1)*h);end for i=1:n2b(i,1)=3*(h2);%初始右端项x(i,1)=0;%初始值endtic%共轭斜量法Ax=br=b-A*x;p=r;for k=1:1000 r1=r; a=(r1*r1)/(A*p)*p); x=x+a*p; r=r-a*A*p; s=(r*r)/(r1*r1); p=r+s*p; if norm(r) u,k,t=CG(9)k = 13t =0.0019u的结果如下表：将u画成三维图 surf(u)将三个程序结果进行比较：对于同一问题，SOR迭代方法迭代次数为33次，迭代时间为0.00219s；块SOR迭代方法迭代次数为25次，迭代时