五点差分格式求解第一边值问题

1、五点差分格式求解第一边值问题一五点差分格式构造的介绍考虑Poisson方程：TOC o 1-5 h z一Au=f(x,y),(x,y)GG,1.1G是xy平面上一有界区域，其边界工为分段光滑曲线，在工满足边值条件：u|t=a(x,y)1.2为解决此问题，做矩形网的差分格式，取定X轴和y轴方向的步长hl和h2,沿x、y方向分别用二阶中心差商代替微商uxx，uyy则得卜u.=ui+ij-2uii+ui-ij+uij+i-2uii+uij-i=f.3hijh12h22ij其中uij表示节点(i，j)上的网函数二、模型问题用五点差分格式求解单位正方形区域的Poisson方程第一边值问题:0 x,y1,

2、d2ud2uAu=+=f.dx2dy2iju|t=,其中T为正方形区域的边界*yj=4u23j=3u1223u22u32j=2u21j=1j=0i=0i=1i=2i=31=图1u7u8u9789u4u5U6uiu2U3图2特别取正方形网格hl=h2=h，利用Yaylor展式，u2u.+ud2u(x.,y.)i+i,jijit,j=rjh2dX2h2d4u(x.,y.)+12dx4已+0(h4)dX4u-2u.+ud2u(x.,y.)h2S4u(x.,y.)i，j+iiji，jt=1J+1J0(h4)TOC o 1-5 h zh22dy212dy4带入上式后，再用出代替琐比,”)，并略去误差项得

3、4u.(u+u+u+u)=h2f.,i,j=1,2，n11.41j1+1,j11,j1,j+11,j11ju“=Pn，u=卩，i=0,1,2，ni,0i,0i,ni,nua.=u=甲，j=0,1,2，n0,j0,jn,jn,j现假定申=0，将1.4式写成矩阵形式即为TU+UT=h2F,l.5n1n1其中U=uij,F=fij,2-11.T=1.eR(nDx(nD-12-三、线性方程组的形成若把正方形的顶点按如图2所示的次序排列，即先按j由小到大，j相同的按i由小到大,这种排列方式叫“自然顺序排列”，用该方式排列得到的线性方程组由下面形状Au=h2f,其中TOC o 1-5 h zi-T+2II

4、in-1n-1n-1-IA=n1I1IT+2In1n1n1A为（N-l）2阶矩阵，如果节点的排列次序改变了，A的形状要改变，但可以证明A的特征值仍不改变。上面的系数矩阵A有这样几个特点：（1）A是块三对角矩阵，共有五条对角线上有非零元素；（2）A是不可约对角占优的；（3）A是正定对称的，而且是稀疏的.容易验证几1的特征值X=2（1COS卩对应的单位特征向量为2jn_sin_nnin少nn2sin(n1)jn)nn利用1.5式，由此即可推出A的特征值为入=入+入=2（2cospncosqn）,1.6pqpqnn对应的特征向量Vpq为矩阵ZpZqT按列“拉直”得到的（n-1）2维向量，即v=pq2

5、(sin_qnzt,sin2qnnnpnzt,.,sin(n】)qnzt)1.7四、方程组的求解由以上知A是不可约对角占优的，利用超松弛迭代逼近Poisson方程第一边值问题的五点差分格式可写成:uk+i=Luk+3(D3L)if,i,j=l,2，n1ij3ijL3=(D3L)1(3R+(13)D),其中3H0是松弛因子,D为A的对角元矩阵-L和-R分别为A的下三角和上三角矩阵，这里松弛因子选为3=3opt=21+J1+p(B)2B为Jacobi迭代矩阵，且B=I-4A所以B矩阵的特征值为怙=14入pq=2吩+吧），p,q=l,2,，n1于是p(B)=cosn=coshnn五、程序设计思想六、

6、附件#include#includeusingnamespacestd;float*one_array_malloc(intn);/一维数组分配float*two_array_malloc(intm,intn);/二维数组分配floatmatrix_category(float*x,intn);intmain()constintMAX=100;最大迭代次数intn,i,j,k;float*a;float*x_0;/初始向量float*x_k;/迭代向量floatprecision;/精度floatw;/松弛因子coutvv输入精度e:;cinprecision;coutvvendlvv输入系数

7、矩阵的阶数，N：;cinn;a=two_array_malloc(n,n+1);coutvvendlvv输入增广矩阵的各值：n;for(i=0;iaij;x_0=one_array_malloc(n);coutvvendlvv输入初始向量:n;for(i=0;ivn;i+)cinx_0i;x_k=one_array_malloc(n);coutvv输入松弛因子w(1vwv2)：n;cinw;floattemp;/迭代过程for(k=0;kMAX;k+)for(i=0;in;i+)temp=0;for(j=0;ji;j+)temp=temp+aij*x_kj;x_ki=ain-temp;temp

8、=0;for(j=i+1;jn;j+)temp=temp+aij*x_0j;x_ki=(x_ki-temp)/aii;x_ki=(1-w)*x_0i+w*x_ki;/求两解向量的差的范数for(i=0;in;i+)x_0i=x_ki-x_0i;if(matrix_category(x_0,n)precision)break;elsefor(i=0;in;i+)x_0i=x_ki;/输出过程if(MAX=k)coutvv迭代不收敛5;coutvv迭代次数为：vvkvvendl;coutvv解向量为：n;for(i=0;ivn;i+)coutxi:x_kiendl;return0;float*one_array_malloc(intn)/一维数组分配float*a;a=(float*)malloc(sizeof(float)*n);returna;float*two_array_malloc(intm,intn)/二维数组分配float*a;inti;a=(float*)malloc(m*sizeof(float*);for(i=0;im;i+)ai