Laplace九点差分格式

中南林业科技大学

本科课程设计说明书

学 院： 理学院

专业年级： 2008级信息与计算科学二班

课 程： 科学计算课程设计

论文题目： Laplace方程九点差分格式

指导教师： 陈宏斌

2011年6月

二维椭圆边值问题的九点差分格式

1问题：Laplace方程

uxx uyy 0, x,y G,

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

u (x,y)(Drichlet边值条件).

在此考虑G为正方形区域，G={(x,y) | a<x<b, a<y<b}.

背景：拉普拉斯方程（Laplace'sequation），又名调和方程、位势方程，是一种偏微分方程。因为由法国数学家拉普拉斯首先提出而得名。求解拉普拉斯方程是电磁学、天文学和流体力学等领域经常遇到的一类重要的数学问题，因为这种方程以势函数的形式描写了电场、引力场和流场等物理对象（一般统称为“保守场”或“有势场”）的性质。

2区域剖分

区域G是一个正方形区域，边界a<x<b,a<y<b. 分别沿x，y轴，在[a,b]上取N+1个节点：a x0 x1 xN b,a y0 y1 yN b. 则步长 h=(b-a)/N;

3九点差分格式

3.1向量形式

利用Taylor展式，有

（3.1）

u(xi 1,yj) 2u(xi,yj) u(xi 1,yj)

2

h x2

68

h4 u(xi,yj)2h6 u(xi,yj)8 o(h),68360 x8! x

（3.2）

2u(xi,yj)

4

h2 u(xi,yj) 12 x4

4

h2 u(xi,yj) 12 y4

u(xi,yj 1) 2u(xi,yj) u(xi,yj 1)

h

4

6

2

6

8

2u(xi,yj) y2

h u(xi,yj)2h u(xi,yj)

o(h8),68360 y8! y

将（3.1），（3.2）两式相加，则得

4u(xi,yj) 4u(xi,yj) h4 6u(xi,yj) 6u(xi,yj) h2 hu(xi,yj) u(xi,yj) 4466 12 x y x y 360 8u(xi,yj) 8u(xi,yj) 2h6 o(h8), 88 8! x y

其中

hu(xi,yj)

u(xi 1,yj) 2u xi,yj u(xi 1,yj)

h

2

u(xi,yj 1) 2u(xi,yj) u(xi,yj 1)

h

2

,

u(xi,yj)

又

2u(xi,yj) x2

2u(xi,yj) y2

0.

6u(xi,yj) x6

6u(xi,yj) y6

2u xi,yj 2u(xi,yj) 6u(xi,yj) 6u(xi,yj) 4 4 224224 x4 y4 x y x y x y 2u xi,yj 2u xi,yj 4 0, 22 22 x y x y

6u(xi,yj) x y

4

2

6u(xi,yj) x2 y4

4u xi,yj 4u xi,yj 2 2u xi,yj 2u xi,yj 4u xi,yj 4u xi,yj 2 2 2 2,442222222 x y x y x y x y x y

而

4u xi,yj x2 y2

u''xx(xi,yj 1) 2u''xx(xi,yj) u''xx(xi,yj 1)

h2

2u xi,yj h4 8u xi,yj h2 4 o(h6)4 262 360 y x12 y x

1

u(xi 1,yj 1) 2u xi,yj 1 u(xi 1,yj 1) 2(u(xi 1,yj) 2u(xi,yj) u(xi 1,yj)) u(xi 1,yj 1) 2u(xi,yj 1) u(xi 1,yj 1) h4

4

4u(xi,yj) 6u(xi,yj 1) h2 6u(xi,yj 1) 6u(xi,yj) 6u(xi,yj 1) h2 6u(xi,yj)1 u(xi,yj 1)

2 2 4446662412 x x x360 x x x 12 x y

8

h4 u(xi,yj) o(h6)26360 x y

1

u(xi 1,yj 1) 2u xi,yj 1 u(xi 1,yj 1) 2(u(xi 1,yj) 2u(xi,yj) u(xi 1,yj)) u(xi 1,yj 1) 2u(xi,yj 1) u(xi 1,yj 1) h4

66888

h2 u(xi,yj)h2 u(xi,yj)1h4 u(xi,yj)h4 u(xi,yj)h4 u(xi,yj) * o(h6)422444622612 x y12 x y1212 x y360 x y360 x y

1

u(xi 1,yj 1) 2u xi,yj 1 u(xi 1,yj 1) 2(u(xi 1,yj) 2u(xi,yj) u(xi 1,yj)) u(xi 1,yj 1) 2u(xi,yj 1) u(xi 1,yj 1) h4

888

1h4 u(xi,yj)h4 u(xi,yj)h4 u(xi,yj) * o(h6)4462261212 x y360 x y360 x y

1

u(xi 1,yj 1) 2u xi,yj 1 u(xi 1,yj 1) 2(u(xi 1,yj) 2u(xi,yj) u(xi 1,yj)) u(xi 1,yj 1) 2u(xi,yj 1) u(xi 1,yj 1) h4

8

h4 u(xi,yj) o(h6)26720 x y

因此

4

8u(xi,yj) 8u(xi,yj) h2 u xi,yj 2h6 8 o(h) hu(xi,yj) u(xi,yj) 8822 8! x y6 x y

1

u(xi 1,yj 1) 2u xi,yj 1 u(xi 1,yj 1) 2(u(xi 1,yj) 2u(xi,yj) u(xi 1,yj)) u(xi 1,yj 1) 2u(xi,yj 1) u(xi 1,yj 1) 6h2

8

8u(xi,yj) 8u(xi,yj) h2h4 u(xi,yj)2h6 8 o(h) u(xi,yj) 8826 8! x y6720 x y

舍去截断误差项，可得Laplace 方程的九点差分格式

1

20ui,j 4ui 1,j 4ui 1,j 4ui,j 1 4ui,j 1 ui 1,j 1 ui 1,j 1 ui 1,j 1 ui 1,j 1 0

6h2

从而可得差分算子 h的截断误差

8

8u(xi,yj) 8u(xi,yj) h2h4 u(xi,yj)2h6 o(h8) Ri,j(u) 88 8! x y6720 x2 y6

40h6

Ri,j(u) M8,

3*8!

其中M8是u的8阶偏导数的绝对值于考虑区域G 1,且u是Laplace 方程的

的光滑解.

3.2 矩阵形式

Au b其中

A1 A 2

A

A2

A1

204 4 204

, A1 ,

A2 4 204

A1 4 20 (N 1)*(N 1) (N 1)*(N 1)

A2

A2

A1A2

A2

4 1

1

4

1 1

41

; 1 4 (N 1)*(N 1)

u u1,u2, ,uN 1 '，

ui ui1,ui2,ui3, ,ui,N 1 ',i 1,2,3, ,N 1;

b b1,b2,b3, ,bN 1 ',

b1 u00 4u01 u02 4u10 u20, u01 4u02 u03, u02 4u03 u04, , u0,N 3 4u0,N 2 u0,N 1, u0,N 2 4u0,N 1 u0,N 4u1,N u2,N ,bi ui 1,0 4ui,0 ui 1,0,0, ,0, ui 1,N 4ui,N ui 1,N ,i 2,3,4, N 2;

bN 1 uN 2,0 4uN 1,0 uN0 4uN1 uN2, uN1 4uN2 uN3, uN2 4uN3 uN4, , uN,N 3 4uN,N 2 uN,N 1, uN,N 2 4uN,N 1 uN 2,N 4uN 1,N uN,N ,

4数值试验

2u 2u x2 y2 0;

u(x,y) ex(siny cosy); 0 x 1,0 y 1.

边值条件：

u(0,y) siny cosy, u(1,y) e(siny cosy), xu(x,0) e,

x u(x,1) e(sin1 cos1).

5编程实现

5.1算法

方程Au b中，A是大型稀疏矩阵矩阵，用块Gauss-Seidel迭代法解；

步骤：1输入区间上限、下限，节点；

2分解矩阵A，A D L U;则迭代矩阵：B (D L) 1U； 3由分块矩阵乘法得到块Gauss-Seidel 迭代法的具体形式：

u(0) 0; (k 1)(k)

A2u2 b1; A1u1

(k 1)k 1)k)

A2ui( A2ui( 11 bi,i 2,3, ,N 2; A1ui

(k 1)(k 1) Au Au1N 12N 2 bN 1;

k 0,1,2, .

4输出.

5.2流程图

5.4编写代码

5.5测试

例1

步长 h=1/6

数值解：

1.3362 1.4653 1.5609 1.6172 1.6393 1.5423 1.6926 1.8067 1.8649 1.8897 1.7216 1.9659 2.1132 2.1590 2.0944 1.9889 2.3380 2.5190 2.5661 2.4093 2.5776 2.8509 3.0488 3.1403 3.1553 精确解：

1.3610 1.5029 1.6031 1.6589 1.6688 1.6078 1.7754 1.8939 1.9598 1.9714 1.8994 2.0974 2.2373 2.3152 2.3290 2.2439 2.4778 2.6431 2.7351 2.7513 2.6508 2.9272 3.1224 3.2312 3.2503 例2

步长h=1/8 数值解：

1.2377 1.3340 1.4182 1.4831 1.5260 1.5498 1.3635 1.4637 1.5598 1.6322 1.6762 1.6973 1.4509 1.6028 1.7269 1.8113 1.8517 1.8508 1.5485 1.7708 1.9349 2.0355 2.0702 2.0351 1.6894 1.9993 2.2064 2.3254 2.3588 2.2917 1.9420 2.3330 2.5632 2.6982 2.7441 2.6806 2.5580 2.7946 3.0092 3.1568 3.2312 3.2330

精确解：

1.2656 1.3783 1.4694 1.5377 1.5819 1.6015 1.4341 1.5618 1.6651 1.7424 1.7926 1.8147 1.6250 1.7697 1.8868 1.9744 2.0313 2.0564 1.8414 2.0054 2.1380 2.2373 2.3017 2.3302 2.0866 2.2724 2.4227 2.5352 2.6082 2.6404 2.3644 2.5749 2.7453 2.8728 2.9555 2.9920 2.6792 2.9178 3.1108 3.2553 3.3490 3.3904

取步长h1=h2=1/10，其数值解和精确解的曲面如下：

1.5621 1.7197 1.8206 1.9311 2.0957 2.4086 3.2172 1.5961 1.8086 2.0494 2.3223 2.6315 2.9819 3.3789

5.6附源程序代码

%Laplace方程九点差分格式 clear;

%a=input('输入区间下限：'); %c=input('输入区间上限：'); N=input('输入节点数：'); a=0;c=1;%N=6; n=10;%迭代次数

%网格剖分 h=(c-a)/N; x=a+h:h:c; y=x;

%解Au=b

%对矩阵b赋值 for p=1:N-1

for q=1:N-1 b(p,q)=0; if p==1

if q==1

b(p,q)=-4*(sin(y(q))+cos(y(q)))-1-(sin(y(2))+cos(y(2)))-4*exp(x(1))-exp(x(2));%边值条件 end

if q==N-1

b(p,q)=-lap9_u(a,y(N-2))-4*lap9_u(a,y(N-1))-lap9_u(a,c)-4*lap9_u(x(1),c)-lap9_u(x(2),c); end

if q>1&q<N-1

b(p,q)=-lap9_u(a,y(q-1))-4*lap9_u(a,y(q))-lap9_u(a,y(q+1)); end end

if p>1&p<N-1 if q==1

b(p,q)=-lap9_u(x(1),a)-4*lap9_u(x(2),a)-lap9_u(x(3),a); end

if q==N-1

b(p,q)=-lap9_u(x(1),c)-4*lap9_u(x(2),c)-lap9_u(x(3),c); end end

if p==N-1 if q==1

b(p,q)=-lap9_u(x(N-2),a)-4*lap9_u(x(N-1),a)-lap9_u(c,a)-4*lap9_u(c,y(1))-lap9_u(c,y(2)); end

if q==N-1

b(p,q)=-lap9_u(x(N-2),c)-4*lap9_u(x(N-1),c)-lap9_u(c,c)-4*lap9_u(c,y(N-1))-lap9_u(c,y(N-2)); end

if q>1&q<N-1

b(p,q)=-lap9_u(c,y(q-1))-4*lap9_u(c,y(q))-lap9_u(c,y(q+1)); end end end end

%对矩阵u赋初值

for p=1:N-1

for q=1:N-1 u0(p,q)=1; end end

%赋值给矩阵A1、A2

b11=-20*ones(1,N-1); b22=4*ones(1,N-2); b33=ones(1,N-2); b44=4*ones(1,N-1);

A1=diag(b11)+diag(b22,1)+diag(b22,-1); A2=diag(b33,1)+diag(b33,-1)+diag(b44);

%用块Gauss-Seidel迭代法解Au=b for m=1:100

u1(1,:)=inv(A1)*(-A2*u0(2,:)'+b(1,:)'); for k=2:N-2

u1(k,:)=inv(A1)*(-A2*u1(k-1,:)'-A2*u0(k+1,:)'+b(k,:)'); end

u1(N-1,:)=inv(A1)*(-A2*u1(N-2,:)'+b(N-1,:)'); u0=u1; end

disp('数值解：'); disp(u1);

%真值

for p=1:N-1

for q=1:N-1

t(p,q)=lap9_u(x(p),y(q)); end end

disp('精确解：'); disp(t);

%画图

x=a+h:h:c-h; y1=x;

[X,Y]=meshgrid(x,y1); subplot(2,1,1); mesh(X,Y,u1); title('数值解曲面');

subplot(2,1,2); surf(X,Y,t);

本文来源：https://www.bwwdw.com/article/m3t1.html