五点菱形格式-拉普拉斯第一边值问题-计算实习报告

计算实习报告姓名：蔡欣麟学号：2013302736一、实习目的1、熟悉偏微分方程数值解的理论知识；2、提高matlab编程能力。3、进一步加深对使用matlab解决相关数学问题的理解。二、实习内容运用matlab实现对Laplace方程第一边值问题的求解。三、算法公式考虑如下区域𝛺̅=𝛺∪𝜕𝛺上求解Poisson方程边值问题Lu=−∆u=f(x,y)∈(a,b)×(c,d)=Ω(1)(∂u∂n+k(x,y)u)𝜕𝛺=𝑟(𝑥,𝑦)的差分方法。用差分方法求解椭圆边值问题时，对区域离散化统一采用矩形网格剖分。对𝛺̅是矩形区域的情况，可将[a,b]等分为N等分，记ℎ1=𝑎𝑏̅̅̅/𝑁；将[c,d]等分为M等分，记ℎ2=𝑐𝑑̅̅̅/𝑀；Ω的离散结果记为Ωℎ，网格节点记为（i,j）,边界结点为i=1和N+1,j=1和M+1的情形，内结点为2≤i≤N，2≤j≤M的情形。在差分方程中，微分方程的离散化处理都采用差商代替微商的方法。为此，要先假设解函数u（x,y）足够光滑，以便利用Taylor展开式做出差商近似。设（i,j）是内节点，为了表现在该结点处的微商，应用Taylor展式，有1ℎ12[𝑢(𝑥𝑖+1,𝑦𝑖)−2𝑢(𝑥𝑖,𝑦𝑖)+𝑢(𝑥𝑖−1,𝑦𝑖)]=𝜕2(𝑥𝑖,𝑦𝑖)𝜕𝑥2+ℎ1212𝜕4𝑢(𝑥𝑖,𝑦𝑖)𝜕𝑥4+ℎ14360𝜕6𝑢(𝑥𝑖,𝑦𝑖)𝜕𝑥6+𝑂(ℎ16)1ℎ22[𝑢(𝑥𝑖,𝑦𝑖+1)−2𝑢(𝑥𝑖,𝑦𝑖)+𝑢(𝑥𝑖,𝑦𝑖−1)]=𝜕2𝑢(𝑥𝑖,𝑦𝑖)𝜕𝑦2+ℎ2212𝜕4𝑢(𝑥𝑖,𝑦𝑖)𝜕𝑦4+ℎ24360𝜕6𝑢(𝑥𝑖,𝑦𝑖)𝜕𝑦6+𝑂(ℎ26)用以上两式中的沿x和y方向的二阶中心差商直接代替方程（1）中的𝑢𝑥𝑥和𝑢𝑦𝑦，就得到−1ℎ12[𝑢𝑖+1,𝑗−2𝑢𝑖𝑗+𝑢𝑖−1,𝑗]−1ℎ22[𝑢𝑖,𝑗+1−2𝑢𝑖𝑗+𝑢𝑖,𝑗−1=𝑓𝑖𝑗(2)由于式（2）中只出现u在点（i,j）及其4个邻点共五个节点处的值；另外，对一切内节点（2）都适用；称（2）为五点差分格式。五点差分格式的截断误差为𝑅𝑖,𝑗=∆𝑢(𝑥𝑖,𝑦𝑖)−∆ℎ𝑢(𝑥𝑖,𝑦𝑖)=𝑂(ℎ2)四、程序设计1、计算实例在区域5.00,5.00:Gyx求解Laplace方程第一边值问题：(要求取x方向与y方向取相同的步长)：yyuxxuxuyuyxu200),5.0(200)5.0,(0)0,(),0(5.0,00(结果：取h＝0.125，利用五点菱形格式，用Gauss-Seidel进行迭代计算得数值结果按先沿x方向后沿y方向为：5.25，12.50，18.75，12.50，25.00，37.50，18.75，37.50，58.25．)2、模块设计(1)网格剖分；(2)近似；(3)Gauss-Seidel迭代；(4)输出结果。3、变量说明q：内点初始值(边界值的算术平均值)u：序号点的函数值U:每次迭代序号点的函数值u1：上次迭代的函数值du：两次迭代函数值的偏差。4、程序清单（附）clear;clc;q=25;u=[0;0;0;0;0;0;q;q;q;25;0;q;q;q;50;0;q;q;q;75;0;25;50;75;100;];fori=1:100forj=1:25U(j,i)=u(j);u1(j)=u(j);endu(7)=(u(2)+u(6)+u(8)+u(12))/4;u(8)=(u(3)+u(7)+u(9)+u(13))/4;u(9)=(u(4)+u(8)+u(10)+u(14))/4;u(12)=(u(7)+u(17)+u(11)+u(13))/4;u(13)=(u(8)+u(18)+u(12)+u(14))/4;u(14)=(u(9)+u(19)+u(13)+u(15))/4;u(17)=(u(12)+u(16)+u(18)+u(22))/4;u(18)=(u(13)+u(17)+u(19)+u(23))/4;u(19)=(u(14)+u(18)+u(20)+u(24))/4;forj=1:25du(j)=u(j)-u1(j);endifnorm(du,1)=0.0000001break;endendfori=1:5forj=1:5u2(i,j)=u(5*(i-1)+j);endendmesh(u2);五、结果分析这样我们就求出了这25个点的精确值，在[0,0.5]×[0,0.5]区域内的所有点都可以用最靠近它的网格点来近似，这样我们就很好地解决了这个偏微分方程问题。六、实习心得通过本学期的计算实习，我对偏微分数值解有了更深的了解。虽然学艺不精，对matlab的掌握还有待提高，但是较以前有了明显进步，应该说是这门课让我在这学期受益良多。此外，这也更坚定了我学好matlab等编程软件以解决与自己专业相关问题的决心。最后，感谢王老师这一学期的陪伴与悉心教导，在这里想真心地对王老师说声，老师您辛苦了！谢谢您!