五点差分格式求解泊松方程的第一边值问题(可编辑)

五点差分格式求解泊松方程的第一边值问题（可编辑）五点差分格式求解泊松方程的第一边值问题

摘要:给出了二维泊松方程在单位正方形上的五点差分格式。并运用线性方程组的古典迭代解法??Jacobi迭代求解出在区域上的数值解。最终绘制数值解的图形。

关键字:泊松方程五点差分格式 Jacobi迭代

有限差分法的介绍

有限差分法是求解偏微分方程的主要数值解法之一;其基本思想是把连续问题离散化,即对求解区域做网格剖分,用有限个网格点代替连续区域;其次将微分算子离散化,从而把微分方程组的问题化为线性方程组的求解问题,解方程组就可以得到原问题在离散点上的近似解。

差分法的步骤:1 对求解域做网格剖分

2 插值函数的选择

3 方程组的建立

4 方程组的求解

五点差分格式的构造

二维泊松方程:

在单位正方形上,在正方形边界上的边界条件.在正方形网格上,就是在上离散化,.对于N3如图1所示:

图1

沿方向分别用二阶中心差商代替

2.1

2.2

1、2式相加可得差分方程:

2.3

利用Taylor展式

可得差分算子的截断误差

其中是方程2.3的光滑解。

由于差分方程2.3中只出现在及其四个邻点上的值见图1的中间的粗的点,所以称为五点差分格式。

由边界条件知道,因而2.3式确定了一组具有个未知量的个线性方程。对应的系数矩阵为对称、不可约对角占优,且对角元为正,因而系数矩阵非奇异,且为对称正定阵。

三、方程组的求解

我们已经知道,利用差分方法解椭圆型方程边值问题归结为解大型线性代数方程组的问题。因为差分格式产生的大型线性代数方程组的系数矩阵中非零元素占的比例小,分布很有规律性。而且通过数值线性代数的学习,我们知道对于大型的稀疏矩阵来说,迭代法是比较好的选择,其程序实现比较简单,迭代过程能自动校正计算过程中的偶然误差,要求计算机的存储相对较少。

本文采用了线性方程组古典迭代解法??Jacobi迭代求解由五点差分格式得到的线性方程组。以下对Jacobi迭代作简要的介绍:

给定3.1

令3.2

其中

3.3

那么3.1可以写成,3.4

其中.若给定初始向量,并代入3.4的右端,就可以计算出一个新的向量,即,

再把代入3.4的右端,又可以得到一个向量;依次类推有,.这就是Jacobi迭代格

式.称为Jacobi迭代的迭代矩阵,称为常数项.

四、算法及流程图

1算法:

输入整数NN可取自2n+1n1,2,3,…构成称数列中的任意数;误差要求e;最大迭代次数M。

输出的近似值;其中

Step1 Set

Step2 For 构造网格点

Set

Step3 For

Set

Step4 For 对网格点按列编号 For ForSet

Step5 For 构造方程组的系数矩阵A和方程组右端项Set b[L]h*h

A[L][L]4If A[L][left] is bounderyThen b[L] +A[L][left]

Else A[L][left]-1

If A[L][right] is bounderyThen b[L] +A[L][right]

Else A[L][right]-1

If A[L][up] is bounderyThen b[L] +A[L][up]

Else A[L][up]-1

If A[L][down] is bounderyThen b[L] +A[L][down]

Else A[L][down]-1

Step6 执行Jacobi迭代输出结果

2流程图:

图4-1

五、数值结果

在给定步长h0.125,将计算机计算出的结果在Matlab中绘制图5-1.如下: 在单位正方形区域上五点差分得到的泊松方程的解图5-1 参考文献:

[1]李荣华冯国忱,微分方程数值解法第三版,高等教育出版社,1995. [2]徐树方高立张平文,数值线性代数,北京大学出版社,2000. [3]James W. Demmel,应用数值线性代数,人民邮电出版社,2007. [4]戴嘉尊邱建贤,微分方程数值解法,东南大学出版社,2002. [5]Arieh lserles ,微分方程数值分析基础教程,清华大学出版社,2005. 六、附录

C程序:

#include

#include

#define N 7

#define h 1/N+1

#define M 100

main

double jacobiint A[N*N][N*N],double b[N*N];

int i,j,L,l,r,up,d,bc0;

double u[N+2][N+2]0; /*包含边界点*/

double v[N*N];

int A[N*N][N*N]0;

double B[N*N];/*初始化矩阵B*/

fori0;iN*N;i++B[i]0.125*0.125; /*按列编号*/ fori1;iN;i++

forj1;jN;j++ forL0;LN*N;L++ v[L]u[i][j]; /*给矩阵A赋值*/

forL0;LN*N;L++ A[L][L]4; ifL/N0 B[L]B[L]+bc; else lN*L/N-1+L%N;

A[L][l]-1; ifL/N+2N B[L]B[L]+bc; else rN*L/N+1+L%N;A[L][r]-1;

ifL%N+2N B[L]B[L]+bc; elseupN*L/N+L%N+1; A[L][up]-1; ifL%N0 B[L]B[L]+bc; else dN*L/N+L%N-1; A[L][d]-1; /*输出矩阵A*/ printf"输

出系数矩阵A: "; printf"\n"; fori0;iN*N;i++

forj0;jN*N;j++printf"%5d ",A[i][j]; printf"\n"; printf"\n";

/*输出矩阵b*/

printf"输出b: "; printf"\n";

fori0;iN*N;i++ printf"%f ",B[i]; printf"\n"; /*调用函数计算结果

*/ v[N*N]jacobiA,B;

system"pause";

double jacobiint A[N*N][N*N],double b[N*N]

double D[N*N][N*N]0;

double L[N*N][N*N]0;

double U[N*N][N*N]0;

double B[N*N][N*N]0;

int i,j,m;

double g[N*N]0;

double x[M][N*N]1;

double v[N*N]; /*给D矩阵赋值*/

fori0;iN*N;i++

/*D[i][i]1/A[i][i];*/D[i][i]0.25; /*输出矩阵D*/printf"输出矩阵D: "; printf"\n"; fori0;iN*N;i++