泊松方程的数值解法-论文.docx

1、泊松方程的数值解法 指导老师：何银年 班级：数学实验班01姓名：宋鹏飞学号：10092014摘要：泊松方程是工程及物理中常见的椭圆型偏微分方程。本文将使用现代数值解法中的有限元法，有限差分法和有限体积法研究泊松方程。关键词：泊松方程 有限元 有限差分 有限体积正文： 泊松方程是数学中一个常见于静电学、机械工程和理论物理的椭圆型偏微分方程，由法国数学家、几何学家及物理学家泊松而得名的。 泊松方程： - 这里是拉普拉斯算子，是欧几里得空间。 一，泊松方程的几种数值解法简介： 泊松方程的数值解有很多种，在这里，我们介绍三种有限元法，有限差分法，有限体积法。有限元法：有限元法实际上就是Ritz-Gal

2、erkin法，下面介绍一下有限元法的具体原理：考虑泊松方程的物理意义，使用最小位能原理，我们知道是位能最小时的的解，位能表达式为： 求的最小值，这让我们想起了线性方程的解法：，要求x的极小值，我们转化成的极值问题，就是对求导，让，最后求出x。同样，对于，我们化成，后求导，让，再使用泛函中的虚功原理，我们就将微分形式变成变分形式： 但是就算如此，无限维的泊松方程还是很难解，我们可以考虑其近似数值，解，使用Ritz-Galerkin法，方法的基本思想是在中选一个N位子空间，然后再上解上述变分问题。总结上述所说的，有限元法可分为有以下几个步骤：（1） 选择单元结构，剖分区域。（2） 构造基函数，形成

3、有限元空间（3） 形成有限元方程，解有限元方程（4） 收敛性及误差估计有限差分法：和有限元法比较，有限差分法就是在第三步中，有限差分将微分形式变成差分形式，其他几步都是一样的。下面给一个差分例子。考虑一个椭圆方程： 方程的差分形式是： 将上述所说的总结一下就推出有限差分法的几个步骤：（1） 对求解区域进行网格剖分（2） 构造逼近微分方程定解问题的差分格式（3） 差分方程的解法（4） 收敛性，稳定性研究有限体积法：其实，有限体积法可以理解为广义差分法。传统的差分法只限于矩形网络，而矩形网络逼近一般误差大，不灵活。所以可以考虑将差分法推广到三角网络，这就是广义差分法，即有限体积法。有线体积法与有限

4、差分法最大区别在于差分法使用了广义导数，并将区域剖分对偶化。二，泊松方程几种数值解法详细介绍：在前面，我们已经给出了有限元法，有限差分法，有限体积法的基本步骤，这里，我们根据这些步骤来对一维和二维的泊松方程进行求解以及误差估计。一， 一维泊松方程的数值解：（1）区域剖分：这里，将之剖分 我们定义 （2）构造有限元空间：为了得到有限元方程和差分方程，我们从Ritz法出发，构造试探函数即基函数，如下：是由组成的函数空间，其中同时，我们定义： （3） 有限元方程，有限差分方程，广义差分方程事实上我们可以通过 ，求来求泊松方程数值解这样我们得到了数值方程如下：有限元方程： 其中有限差分方程：有限体积法

5、方程：（4）解方程组：根据这三个方程就可以得到一组线性代数方程，从而可得到三个矩阵 其中 有限元法 有限差分法 有限体积法 解此矩阵即可可得到：（4）稳定性及误差估计：对于稳定性，可以使用误差估计，可以使用二，二维泊松方程的数值解（1）区域剖分：对于x轴，剖分和一维的一样，对于y轴，我们定义：并设定对于每个子区域，把它分为六个小区域 如图：（2）构造有限元空间：是由构成的函数空间，其中定义如下边界情况一样定义（3）有限元方程，有限差分方程，广义差分方程：事实上我们可以通过求来求泊松方程数值解这样我们得到了数值方程如下：有限元方程：有限差分方程： 其中：有限体积法方程： 其中：（4）解方程组：根

6、据这三个方程就可以得到一组线性代数方程，从而可得到三个矩阵其中有限元法；有限差分法；有限体积法 解此矩阵即可。解的结果通式如下：对于有限元和有限体积区域剖分法是一样的所以，A相同： 其中：有限体积法的区域剖分和前两种方法不同，A也不同：其中：（5）稳定性及误差估计：稳定性可以使用：其中误差估计使用：三，数值实验例子：例：使用差分法解： 其中：算法结果：U = 0.009396957688515 0.018645427876649 0.027597294205442 0.036105181671115 0.044022824985826 0.042826878058501 0.030908304

7、319090 0.016662475915615最大的误差是0.000015如图：误差图：参考文献：1,偏微分方程现代数值解法，马逸尘，梅立泉，王阿霞编著，2006年10月，科学出版社2,偏微分方程数值解法，李荣华，2005年5月3,数值分析，李乃成，梅立泉，2011年9月，科学出版社4，精通matlab，葛哲学，电子工业出版社附录：数值matlab程序：syms x;a=0;b=1;%区间界点N=20;%将区间划分的等分，这里控制！h=(b-a)/N; %这里确定步长f=sin(x)/sin(1);%这是f函数 g=(h2)*f;value_of_g=zeros(N-1,1);%这是（h2)

8、*fdiag_0=zeros(N-1,1);%确定A的对角元diag_1=zeros(N-2,1);%确定A的偏离对角的上对角元diag_2=zeros(N-2,1);%确定A的偏离对角的下对角元X=a:h:b;y_a=0;%边界条件y_b=0;%边界条件for j=2:N diag_0(j-1)=2;end %获取对角元素 for j=3:N diag_1(j-2)=-1;end %获取A的第三条对角for j=2:N-1 diag_2(j-1)=-1;end %获取A的第二条对角 for j=2:N value_of_g(j-1)=subs(g,x,X(j);end %获取G值A=diag

9、(diag_0)+diag(diag_1,1)+diag(diag_2,-1);%组装矩阵Aformat longU=Avalue_of_g %差分解%fprintf(%11.5f,U)fprintf(n);dx=X(2:N);precise_value=sin(dx)./sin(1)-dx %精确解，为已知%fprintf(%11.5f,precise_value)deta=U-precise_value; %误差deta_max=max(abs(deta); %最大误差fprintf(最大的误差是%fn,deta_max)plot(X(2:N),U,b*,X(2:N),precise_value,r-)%差分解与精确解对比图gridlegend(差分解,精确解,Location,NorthEast)xlabel(x)；ylabel(y)figure();plot(X(2:N),deta) %误差图grid；xlabel(x)；ylabel(deta)