数值求解如下正方形域上的Poisson方程边值问题
2u2u22f(x,y)1,0x,y1yx u(0,y)u(1,y)y(1y),0y10x1u(x,0)u(x,1)0,
二、用椭圆型第一边值问题的五点差分格式得到线性方程组为
4ui,jui1,jui1,jui,j1ui,j1h2fij
1i,jNu0,j?,uN1,j?,ui,0?,ui,N1?,0i,jN1
写成矩阵形式Au=f。其中 A11IAIA22I
IANNv1vu2vNb1bf2bN4114Aii114v1(u1,1,u2,1,...,uN,1)T,v2(u1,2,u2,2,...,uN,2)T,......,vN(u1,N,u2,N,...,uN,N)Tb1h2(f1,1,f2,1,...,fN,1)T?,b2h2(f1,2,f2,2,...,fN,2)T?,......,bNh2(f1,N,f2,N,...,fN,N)T?1h,取N9或99,N1fi,j1,i,j1,2,...,N
三 、编写求解线性方程组Au=f的算法程序,用下列方法编程计算,并比较计算速度。 1. 2. 3. 4. 5. 6.
用Jacobi迭代法求解线性方程组Au=f。 用块Jacobi迭代法求解线性方程组Au=f。
用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。 用块SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。 用Gauss-Seidel迭代法求解线性方程组Au=f。 用块Gauss-Seidel迭代法求解线性方程组Au=f。
则h0.1或0.01四、上机报告要求
1.简述方法的基本原理。 2.程序中要加注释。
3.对程序中的主要变量给出说明。 4.附原程序及计算结果。
5.迭代法效率分析:比较迭代法的效率,记录所需的迭代次数,所花费的计算时间。 6.此题可分组做,每组1-3人。需要注明每人完成的任务。
附参考程序:
1、用Gauss-seidel迭代法求解线性方程组Au=f,具体程序如下: function [u,k]=xsgs(n)
f(2:n+1,2:n+1)=(n+1)^(-2)*2; u=zeros(n+2,n+2); e=0.000000001;
for k=1:1000 %迭代求解 er=0;
for j=2:n+1 for i=2:n+1 Ub=u(i,j);
u(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-u(i,j)); %估计误差 end end
if er/n^2 % xsbgs:用块Gauss-seidel迭代法求解线性方程组A*u=f % u: 方程组的解; k: 迭代次数; n: 非边界点数 % a: 方程组系数矩阵Aii的下对角线元素; b: 方程组系数矩阵Aii的主对角线元素; % c: 方程组系数矩阵Aii的上对角线元素;d: 追赶法所求方程的右端向量; % e: 允许误差界; er:迭代误差; f=2*1/(n+1)^2*ones(n+2,n+2); a=-1*ones(1,n); b=4*ones(1,n);c=-1*ones(1,n); u=zeros(n+2,n+2);e=0.00001; 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) ; x=zg(a,b,c,d); % 用追赶法求解 u(2:n+1,j)=x'; er=er+norm(Ub-u(:,j),1); end if er/n^2 >> n=9; b(2:n+1,2:n+1)=0.02; U=zeros(n+2,n+2); e=0.000000001; tic; for k=1:1000 %迭代求解 er=0; for j=2:n+1 for i=2:n+1 Ub=U(i,j); U(i,j)=(U(i-1,j)+U(i+1,j)+U(i,j-1)+U(i,j+1)+b(i,j))/4; er=er+abs(Ub-U(i,j)); %估计误差 end end if er/n^2 因篇幅问题不能全部显示,请点此查看更多更全内容