Laplace九点差分格式

中南林业科技大学本科课程设计说明书学院：理学院专业年级：2008级信息与计算科学二班课程：科学计算课程设计论文题目：Laplace方程九点差分格式指导教师：陈宏斌2011年6月Laplace方程九点差分格式共10页第1页二维椭圆边值问题的九点差分格式1问题：Laplace方程,,,0GyxuuyyxxG是xy平面上一有界区域，其边界为分段光滑曲线.在上u满足下列边值条件：)边值条件Drichlet)(,(yxu.在此考虑G为正方形区域，G={(x,y)|axb,ayb}.背景：拉普拉斯方程（Laplace'sequation），又名调和方程、位势方程，是一种偏微分方程。因为由法国数学家拉普拉斯首先提出而得名。求解拉普拉斯方程是电磁学、天文学和流体力学等领域经常遇到的一类重要的数学问题，因为这种方程以势函数的形式描写了电场、引力场和流场等物理对象（一般统称为“保守场”或“有势场”）的性质。2区域剖分区域G是一个正方形区域，边界axb,ayb.分别沿x，y轴，在[a,b]上取N+1个节点：byyyabxxxaNN1010,.则步长h=(b-a)/N;3九点差分格式3.1向量形式利用Taylor展式，有（3.1）),(),(!82),(360),(12),(),(),(2),(888666444222211hoxyxuhxyxuhxyxuhxyxuhyxuyxuyxujijijijijijiji（3.2）Laplace方程九点差分格式共10页第2页),(),(!82),(360),(12),(),(),(2),(888666444222211hoyyxuhyyxuhyyxuhyyxuhyxuyxuyxujijijijijijiji将（3.1），（3.2）两式相加，则得),(),(),(!82),(),(360),(),(12),(),(8888866666444442hoyyxuxyxuhyyxuxyxuhyyxuxyxuhyxuyxujijijijijijijijih其中,),(),(2),(),(,2),(),(211211hyxuyxuyxuhyxuyxuyxuyxujijijijijijijih.0),(),(),(2222yyxuxyxuyxujijiji又426246222244446666),(),(),(,),(),(yxyxuyxyxuyyxuxyxuyxyyxuxyxujijijijijiji,0,,),(),(2222224426246yyxuxyxuyxyxyxuyxyxujijijiji,,2,2,,,,224224222222224444yxyxuyxyxuyyxuxyxuyxyyxuxyxujijijijijiji而)(,360,12),(''),(''2),('',6268422442211224hoxyyxuhxyxuyhhyxuyxuyxuyxyxujijijixxjixxjixxjiLaplace方程九点差分格式共10页第3页)(),(360),(12),(),(2),(360),(),(2),(121),(),(2),()),(),(2),((2),(,2),(1662844262616666162416444141111111111114hoyxyxuhyxyxuhxyxuxyxuxyxuhxyxuxyxuxyxuyxuyxuyxuyxuyxuyxuyxuyxuyxuhjijijijijijijijijijijijijijijijiji)(),(360),(360),(12*121),(12),(12),(),(2),()),(),(2),((2),(,2),(16628426844484426224621111111111114hoyxyxuhyxyxuhyxyxuhyxyxuhyxyxuhyxuyxuyxuyxuyxuyxuyxuyxuyxuhjijijijijijijijijijijijijiji)(),(360),(360),(12*121),(),(2),()),(),(2),((2),(,2),(166284268444841111111111114hoyxyxuhyxyxuhyxyxuhyxuyxuyxuyxuyxuyxuyxuyxuyxuhjijijijijijijijijijijiji)(),(720),(),(2),()),(),(2),((2),(,2),(1662841111111111114hoyxyxuhyxuyxuyxuyxuyxuyxuyxuyxuyxuhjijijijijijijijijiji因此),(),(jijihyxuyxu2242,6yxyxuhji)(),(),(!82888886hoyyxuxyxuhjiji),(),(2),()),(),(2),((2),(,2),(611111111111112jijijijijijijijijiyxuyxuyxuyxuyxuyxuyxuyxuyxuh62842),(7206),(yxyxuhhyxujiji)(),(),(!82888886hoyyxuxyxuhjiji舍去截断误差项，可得Laplace方程的九点差分格式0444420611,11,11,11,11,1,,1,1,2jijijijijijijijijiuuuuuuuuuh从而可得差分算子h的截断误差62842,),(7206)(yxyxuhhuRjiji)(),(),(!82888886hoyyxuxyxuhjiji,!8*340)(86,MhuRji其中M8是u的8阶偏导数的绝对值于考虑区域G的上确界，,1且u是Laplace方程的Laplace方程九点差分格式共10页第4页的光滑解.3.2矩阵形式;bAu其中,)1)*(1(NN1221221221AAAAAAAAAAA,20442044204420)1)*(1(NN1A;4114114114)1)*(1(2NNA'1N21u,,u,uu，;1,,3,2,1,',,,,1,321NiuuuuNiiiiiu'1N321b,,b,b,bb,,44,4,,4,4,44,2,1,01,02,01,02,03,004030203020120100201001NNNNNNNNuuuuuuuuuuuuuuuuuuub;2,4,3,2,4,0,,0,4,1,,10,10,0,1NiuuuuuuNiNiNiiiiib,44,4,,4,4,44,,1,21,2,1,2,3,4323212100,10,21NNNNNNNNNNNNNNNNNNNNNNNNNNNNuuuuuuuuuuuuuuuuuuub4数值试验Laplace方程九点差分格式共10页第5页.10,10);cos(sin),(;02222yxyyeyxuyuxux边值条件：).1cos1(sin)1,(,)0,(),cos(sin),1(,cossin),0(xxexuexuyyeyuyyyu5编程实现5.1算法方程bAu中，A是大型稀疏矩阵矩阵，用块Gauss-Seidel迭代法解；步骤：1输入区间上限、下限，节点；2分解矩阵A，;ULDA则迭代矩阵：UL)(DB1；3由分块矩阵乘法得到块Gauss-Seidel迭代法的具体形式：.,2,1,0;;2,,3,2,;;01)1(22)1(11)(12)1(12)1(11)(22)1(11)0(kbuAuANibuAuAuAbuAuAuNkNkNikikikikk4输出.5.2流程图Laplace方程九点差分格式共10页第6页开始输入区间，节点数N给出迭代次数n对向量b赋值，对向量u赋初值u0；对矩阵A1，A2赋值m=1mn+1k=2kN-1YA1*u1(N-1)=-A2*u1(N-2)+b(N-1)A1*u1(1)=-A2*u2(2)+b(1);k=k+1A1*u1(k)=-A2*u1(k-1)-A2*u0(k+1);K=k+1YNm=m+1输出数值解u1N结束5.4编写代码5.5测试例1步长h=1/6Laplace方程九点差分格式共10页第7页数值解：1.33621.46531.56091.61721.63931.54231.69261.80671.86491.88971.72161.96592.11322.15902.09441.98892.33802.51902.56612.40932.57762.85093.04883.14033.1553精确解：1.36101.50291.60311.65891.66881.60781.77541.89391.95981.97141.89942.09742.23732.31522.32902.24392.47782.64312.73512.75132.65082.92723.12243.23123.2503例2步长h=1/8数值解：1.23771.33401.41821.48311.52601.54981.56211.36351.46371.55981.63221.67621.69731.71971.45091.60281.72691.81131.85171.85081.82061.54851.77081.93492.03552.07022.03511.93111.68941.99932.20642.32542.35882.29172.09571.94202.33302.56322.69822.74412.68062.40862.55802.79463.00923.15683.23123.23303.2172精确解：1.26561.37831.46941.53771.58191.60151.59611.43411.56181.66511.74241.79261.81471.80861.62501.76971.88681.97442.03132.05642.04941.84142.00542.13802.23732.30172.33022.32232.08662.27242.4227