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

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

一五点差分格式构造的介绍

考虑Poisson方程：

u=f(x,y)， (x,y)∈G, 1.1

G是xy平面上一有界区域，其边界τ为分段光滑曲线，在τ满足边值条件：

u|τ=α(x,y) 1.2

为解决此问题，做矩形网的差分格式，取定x轴和y轴方向的步长h1和h2，沿x 、y方向分别用二阶中心差商代替微商u xx，uyy 则得

huij= [ui+1,j 2uij+ui 1,jui,j+1 2uij+ui,j 1

h12h22fij 1.3

其中uij表示节点 (i,j) 上的网函数

由于差分格式1.3只出现u在(i,j)及其四个邻点上的值，称为五点差分格式，如图示：

二、模型问题

用五点差分格式求解单位正方形区域的Poisson方程第一边值问题：

2u 2u u=+=fij , 0< , <1 , u|τ=φ ,

其中τ为正方形区域的边界

图 1 图 2

特别取正方形网格h1=h2=h,利用Yaylor展式，

ui+1,j 2uij+ui 1,j 2u(xi,yj) h2 4u(xi,yj)4=++O(h) ui,j+1 2uij+ui,j 1 2u(xi,yj) h2 4u(xi,yj)4=+O(h) 带入上式后，再用uij代替u(xi,yj), 并略去误差项得

4uij ui+1,j+ui 1,j+ui,j+1+ui,j 1 =h2fij ,i,j=1,2,…,n 1 1.4

ui,0=φi,0, ui,n=φi,n, i=0,1,2…,n

u0,j=φ0,j, un,j=φn,j, j=0,1,2…,n

现假定φ=0，将1.4式写成矩阵形式即为

Tn 1U+UTn 1=h2F, 1.5 其中U=[uij], F=[fij],

2 1 1 Tn 1= ∈R n 1 × n 1 1 12

三、线性方程组的形成

若把正方形的顶点按如图2所示的次序排列，即先按j由小到大，j相同的按i 由小到大，这种排列方式叫“自然顺序排列”，用该方式排列得到的线性方程组由下面形状

Au=h2f ,

其中

Tn 1+2In 1 In 1 In 1 A=

In 1 In 1Tn 1+2In 1

A为(N-1)2阶矩阵，如果节点的排列次序改变了，A的形状要改变，但可以证明A 的特征值仍不改变。

上面的系数矩阵A有这样几个特点：

（1） A是块三对角矩阵，共有五条对角线上有非零元素；

（2） A是不可约对角占优的；

（3） A是正定对称的，而且是稀疏的.

容易验证Tn-1的特征值λj=2(1 cos),对应的单位特征向量为 n

Tjπ2jπ22jπ2(n 1)jπzj= sin, sin,…, sin 利用1.5式，由此即可推出A的特征值为

λpq=λp+λq=2(2 cospπn cosqπn, 1.6

对应的特征向量vpq为矩阵zpzqT按列“拉直”得到的(n-1)2 维向量，即

vpq= sinn2 qπnzp,sinT 2qπnzp,…,sinT (n 1)qπnzp 1.7 T

四、方程组的求解

由以上知A是不可约对角占优的， 利用超松弛迭代逼近Poisson方程第一边值问题的五点差

分格式可写成:

1uk+1=Lωukij+ω D ωL f，i,j=1,2…,n 1 ij

Lω= D ωL 1(ωR+(1 ω)D),

其中ω≠0是松弛因子，D为A的对角元矩阵 L 和 R分别为A的下三角和上三角矩阵 ，

=ωopt=这里松弛因子选为ωB为Jacobi迭代矩阵，且B=I A， 1所以B矩阵的特征值为

μpq=1 11pπ

4λpq=2(cosn+cosqπ

n) ,

于是ρ B =cosπ

ncoshπ

五、程序设计思想

4p,q=1,2,…,n 1

六、附件

#include <iostream>

#include <cmath>

using namespace std;

float *one_array_malloc(int n); //一维数组分配 float **two_array_malloc(int m,int n); //二维数组分配 float matrix_category(float* x,int n);

int main()

{

const int MAX=100;//最大迭代次数

int n,i,j,k;

float** a;

float* x_0; //初始向量

float* x_k; //迭代向量

float precision; //精度

float w; //松弛因子

cout<<"输入精度e:";

cin>>precision;

cout<<endl<<"输入系数矩阵的阶数，N："; cin>>n;

a=two_array_malloc(n,n+1);

cout<<endl<<"输入增广矩阵的各值：\n";

for(i=0;i<n;i++)

{

for(j=0;j<n+1;j++)

{

cin>>a[i][j];

}

}

x_0=one_array_malloc(n);

cout<<endl<<"输入初始向量:\n";

for(i=0;i<n;i++)

{

cin>>x_0[i];

}

x_k=one_array_malloc(n);

cout<<"输入松弛因子w (1<w<2)：\n";

cin>>w;

float temp;

//迭代过程

for(k=0;k<MAX;k++)

{

for(i=0;i<n;i++)

{

temp=0;

for(j=0;j<i;j++)

{

temp=temp+a[i][j]*x_k[j];

}

x_k[i]=a[i][n]-temp;

temp=0;

for(j=i+1;j<n;j++)

{

temp=temp+a[i][j]*x_0[j];

}

x_k[i]=(x_k[i]-temp)/a[i][i];

x_k[i]=(1-w)*x_0[i]+w*x_k[i];

}

//求两解向量的差的范数

for(i=0;i<n;i++)

{

x_0[i]=x_k[i]-x_0[i];

}

if(matrix_category(x_0,n)<precision)

{

break;

}

else

{

for(i=0;i<n;i++)

{

x_0[i]=x_k[i];

}

}

}

//输出过程

if(MAX==k)

{

cout<<"迭代不收敛\n";

}

cout<<"迭代次数为："<<k<<endl;

cout<<"解向量为：\n";

for(i=0;i<n;i++)

{

cout<<"x"<<i<<": "<<x_k[i]<<endl;

}

return 0;

}

float *one_array_malloc(int n) //一维数组分配 {

float *a;

a=(float *)malloc(sizeof(float)*n);

return a;

}

float **two_array_malloc(int m,int n) //二维数组分配 {

float **a;

int i;

a=(float **)malloc(m*sizeof(float *));

for (i=0;i<m;i++)

{

a[i]=(float *)malloc(n*sizeof(float)); }

return a;

}

float matrix_category(float* x,int n)

{

int i;

float temp=0;

for(i=0;i<n;i++)

{

temp=temp+fabs(x[i]);

}

return temp;

}