Laplace九点差分格式

1、laplace九点差分格式 中南林业科技高校 本科课程设计说明书 学 院： 理学院 专业班级： 2021级信息与计算科学二班 课 程： 科学计算课程设计 论文题目： laplace方程九点差分格式 指导老师： 陈宏斌 2021年6月 二维椭圆边值问题的九点差分格式 1问题：laplace方程 uxx uyy 0, x,y g, g是xy平面上一有界区域，其边界 为分段光滑曲线. 在 上u满意下列边值条件： u (x,y)(drichlet边值条件). 在此考虑g为正方形区域，g=(x,y) | axb, ayb. 背景：拉普拉斯方程（laplace'sequation），又名调和方程、

2、位势方程，是一种偏微分方程。由于由法国数学家拉普拉斯首先提出而得名。求解拉普拉斯方程是电磁学、天文学和流体力学等领域常常遇到的一类重要的数学问题，由于这种方程以势函数的形式描写了电场、引力场和流场等物理对象（一般统称为“保守场”或“有势场”）的性质。 2区域剖分 区域g是一个正方形区域，边界axb,ayb. 分别沿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(

3、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

4、,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,y

5、j) 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

6、(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

7、(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

8、(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

9、 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) 8

10、u(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

11、, ,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

12、 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由分块矩阵乘

13、法得到块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.98

14、89 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.5

15、498 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.60

16、15 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.

17、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(

18、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 q1qn-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 p1pn-1 if q=1 b(p,q)

19、=-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 q1qn-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 %赋值给