下载文档原格式
泊松方程的有限差分法的MATLAB实现作者：冯立伟徐涛屈福志来源：《电脑知识与技术》2017年第13期摘要：泊松方程是物理及工程应用领域中一类非常重要的方程，研究其数值求解方法具有重要意义。
给出了使用有限差分法求解泊松方程的计算方法，并讨论了使用MATLAB编写计算程序，使用数值算例和静电场实例进行了数值实验，实验结果与理论一致，检验了算法的有效性。
关键词：泊松方程；五点差分格式；有限差分法中图分类号：TP311 文献标识码：A 文章编号：1009-3044（2017）13-0233-031概述物理过程，都可用椭圆型方程来描述。
其中最典型的方程是泊松（Poisson）方程。
传热学中带有稳定热源或内部无热源的稳定温度场的温度分布、流体动力学中不可压缩流体的稳定无旋流动、弹性力学中平衡问题及电磁学中静电场的电势等均满足泊松方程，泊松方程也是数值网格生成技术所遵循的基本方程。
因此，研究其数值求解方法具有重要意义。
MATLAB是目前应用最广泛的科学和工程计算软件。
MATLAB基于矩阵运算，具有强大数值运算能力，是方便实用、功能强大的数学软件；同时，MATLAB具有强大的图形绘制功能，用户只需提供绘图数据和指定绘图方式，用很少的程序指令就可得到将计算结果转化为直观、形象的图像。
使用MAT-LAB求解微分方程已有大量的研究。
因此，近些年来，越来越多的人开始使用MATLAB来求解泊松方程。
利用MAT-LAB强大的数值计算能力和图形绘制技术，可以实现使用差分法求解泊松方程并绘制出数值解的二维、三维图像，从而可以更好地理解泊松方程解的物理意义。
本文讨论使用差分法通过MATLAB编程求解二维矩形区域上的泊松方程，并使用两个算例进行检验和对结果进行分析。
边界条件为将未知解函数在内部节点上的值按行排列，组成解向量为：3差分格式的求解为了便于使用MATLAB编写程序，将差分方程转化为矩阵形式：4数值实验算例1：为了分析和比较差分格式在不同步长下的结果，使用2范数意义下的绝对误差和相对误差作为评价指标，表1给出了步长h=0.01取不同值的绝对误差和相对误差从表1可看出随着网格步长h的减小数值解的绝对误差和相对误差在变小。
高数泊松方程全文共四篇示例，供读者参考第一篇示例：高等数学中的泊松方程是一个非常重要的概念，它在物理学、工程学等领域都有着广泛的应用。
泊松方程是一个偏微分方程，它描述了一个标量函数的拉普拉斯算子（拉普拉斯算子是二阶空间导数的和）等于一个给定函数的情况。
泊松方程在物理学中的应用特别广泛，比如在电磁学中描述电势分布、在热传导中描述温度分布等等。
我们来看一下泊松方程的数学表达式：设函数u(x,y,z)在空间区域D内有定义，记\frac{\partial^2 u}{\partial x^2}、\frac{\partial^2 u}{\partial y^2} 、\frac{\partial^2 u}{\partial z^2} 分别为u对x、y、z的二阶偏导数，那么泊松方程可以表示为：其中f(x,y,z)是给定的函数。
可以看出，泊松方程是一个二阶偏微分方程，描述了一个标量函数的拉普拉斯算子等于一个给定函数的情况。
泊松方程在物理学中的应用非常广泛，其中最为著名的应用就是在电磁学中描述电势分布。
根据麦克斯韦方程组，电荷在空间中分布会产生电场，而电场又会引起电荷的移动，形成电流。
泊松方程可以描述电势分布，即描述电场的强度和方向。
通过求解泊松方程，我们可以得到电势分布的解析表达式，从而进一步推导出电场分布、电流分布等物理量。
另一个比较常见的应用就是在热传导中描述温度分布。
热传导是物体内部的热量传递过程，它遵循热传导方程。
如果我们知道物体表面的温度分布，可以通过泊松方程求解得到物体内部的温度分布。
这对于工程设计和热力分析非常重要，可以帮助我们优化散热结构，提高能效等。
除了电磁学和热传导，泊松方程还有很多其他的应用，比如在流体力学中描述流场的速度分布、在弹性力学中描述物体变形的表面位移等等。
泊松方程是一个非常重要且有着广泛应用的数学工具，它帮助我们理解和描述自然界中许多复杂的现象。
在实际应用中，泊松方程的求解并不是一件容易的事情。
泊松方程泊松方程只得是数学中一个常见于静电学、机械工程和理论物理的偏微分方程，因法国数学家、几何学家及物理学家泊松而得名的。
在数学以及物理中，拉普拉斯算子或是拉普拉斯算符（Laplace operator 或Laplacian）是一个微分算子，通常写成Δ或；这是为了纪念皮埃尔-西蒙·拉普拉斯而命名的。
拉普拉斯算子有许多用途，此外也是椭圆型算子中的一个重要例子。
在物理中，常用于波方程的数学模型、热传导方程以及亥姆霍兹方程。
在静电学中，拉普拉斯方程和泊松方程的应用随处可见。
在量子力学中，其代表薛定谔方程式中的动能项。
在数学中，经拉普拉斯算子运算为零的函数称为调和函数；拉普拉斯算子是霍奇理论的核心，并且是德拉姆上同调的结果。
泊松方程成立的条件泊松首先在无引力源的情况下得到泊松方程,△Φ=0（即拉普拉斯方程）；当考虑引力场时,有▽Φ=f（f为引力场的质量分布）.后推广至电场磁场,以及热场分布.该方程通常用格林函数法求解,也可以分离变量法,特征线法求解.泊松方程的物理内涵泊松方程可以看做是不可压缩的流体运动方程。
方程的意义相当于穿过任意封闭曲面的液体的流量等于曲面内所包含的流体源产生液体的总量。
对于电动力学中静电场，电场强度相当于流密度，净电荷相当于流体源电动力学中电场对空间坐标的二次导数与空间内电荷量成正比。
半导体中的泊松方程泊松方程表明电荷产生电场：电位的二阶导数与电荷密度成正比。
近似条件：PIN结中无载流子即全部耗尽，施主和受主完全电离。
PIN结的泊松方程：(0<x<Xn)d^2V(x)/dx^2=-Nd/ε，(-Xp<x<0)d^2V(x)/dx^2=-Na/ε边界条件E(0)=E(Xn)=-dV(x)/dx(x=-Xp,Xn)=0，V(x=-Xp)=0,V(x=Xn)=0 将上面的式子一次积分（注意符号）带入边界条件就能得出电场的分布，再次积分就能得出电势的分布。
椭圆方程的公式椭圆方程是数学中一个非常重要的概念，它在物理、工程、计算机科学等领域都有广泛的应用。
本文将介绍椭圆方程的公式及其应用。
一、椭圆方程的定义椭圆方程是一个二元二次方程，其一般形式为：Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0其中A、B、C、D、E、F均为实数，且A、C不同时为0。
二、椭圆方程的标准形式椭圆方程可以通过变量替换和平移来化为标准形式：(x-x0)^2/a^2 + (y-y0)^2/b^2 = 1其中(x0,y0)为椭圆中心点坐标，a、b为椭圆长轴和短轴的长度。
三、椭圆方程的参数椭圆方程的参数包括中心坐标、长轴和短轴长度、离心率等。
1. 中心坐标：椭圆的中心坐标为(x0,y0)。
2. 长轴和短轴长度：长轴的长度为2a，短轴的长度为2b。
3. 离心率：椭圆的离心率为e，e的值介于0和1之间，表示椭圆长轴与短轴长度之比。
四、椭圆方程的性质1. 对称性：椭圆方程具有关于x轴和y轴的对称性。
2. 焦点和直径：椭圆方程有两个焦点F1和F2，它们之间的距离为2c，c^2=a^2-b^2。
椭圆的长轴是过焦点F1和F2的直径。
3. 弦和法线：椭圆方程上任意一点P的切线与椭圆长轴的夹角是β，法线与椭圆长轴的夹角是α。
弦是连接椭圆上任意两点的线段，弦的中垂线与长轴的夹角是β/2，法线与弦的夹角是α-β/2。
五、椭圆方程的公式1. 椭圆方程的离心率公式：e=sqrt(1-b^2/a^2)2. 椭圆焦点的坐标公式：F1(x0-c,y0)，F2(x0+c,y0)3. 椭圆长轴和短轴长度公式：a^2=c^2+b^2b^2=a^2-c^24. 椭圆周长公式：C=4aE(e)其中E(e)是第二类椭圆积分，可以用级数或逼近公式计算。
5. 椭圆面积公式：S=πab六、椭圆方程的应用椭圆方程在物理、工程、计算机科学等领域都有广泛的应用，以下是一些例子：1. 圆轨道的近似：当椭圆的离心率e足够小时，它近似为一个圆，因此可以用椭圆方程来描述圆形轨道。
《微分方程数值解》大作业（一）——椭圆型方程编程计算：采用五点差分格式求如下椭圆型方程2222uu x y f (x,y),(x,y);∂∂∂∂--=∈Ω其中f (x,y)、Ω及边条件为：1． f (x,y)0,= (1,2)(0,1)Ω=⨯， 且边条件如下：222u(x,0)2ln x,u(x,1)ln(x 1)1x 2;u(1,y)ln(1y ),u(2,y)ln(4y ),0y 1.⎧==+<<⎪⎨=+=+<<⎪⎩， 问题存在精确解为： 22(,)ln()u x y x y =+2．f (x,y)4,=- (0,1)(0,2)Ω=⨯，且边条件如下：2222u(x,0)x ,u(x,2)(x 2)0x 1;u(0,y)y ,u(1,y)(y 1),0y 2.⎧==-<<⎪⎨==-<<⎪⎩， 问题存在精确解为： 2(,)()u x y x y =-3．f (x,y)cos(x y)cos(x y),=++- (0,)(0,)2πΩ=π⨯，且边条件如下：u(x,0)cos x,u(x,)00x ;2u(0,y)cos y,u(,y)cos y,0y .2π⎧==<<π⎪⎪⎨π⎪=π=-<<⎪⎩， 问题存在精确解为： (,)cos cos u x y x y =.代码：主函数1，差分解function g=fivepoints(x1,x2,y1,y2,M,N)%变步长法h=(x2-x1)/M; %横轴步长k=(y2-y21/N; %纵轴步长m=M-1;n=N-1;h1=h^2;r=h1/k^2; %五点中的上下两个点的系数t=2+2*r; %五点中的中心点的系数x=x1+(x2-x1)*(0:M)/M; %x，y向量表示横纵坐标y=y1+(y2-y1)*(0:N)/N;a=zeros(m*n,m*n);b=zeros(m*n,1);%初始化a，b矩阵，a为系数矩阵%内部的（m-2）*（n-2）个点for i=2:m-1for j=2:n-1a(i+(j-1)*m,:)=[zeros(1,i-1+(j-2)*m) -r zeros(1,m-2) -1 t -1 zeros(1,m-2) -r zeros(1,(n-j)*m-i)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1));endend%下边缘j=1;for i=2:m-1a(i+(j-1)*m,:)=[zeros(1,i-2) -1 t -1 zeros(1,m-2) -r zeros(1,(n-j)*m-i)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*bottom(x(i+1));end;%右边缘i=m;for j=2:n-1a(i+(j-1)*m,:)=[zeros(1,(j-1)*m-1) -r zeros(1,m-2) -1 t zeros(1,m-1) -r zeros(1,(n-j)*m-i)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+right(y(j+1));end%上边缘j=n;for i=2:m-1a(i+(j-1)*m,:)=[zeros(1,i-1+(j-2)*m) -r zeros(1,m-2) -1 t -1 zeros(1,m-i-1)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*top(x(i+1));end%左边缘i=1;for j=2:n-1a(i+(j-1)*m,:)=[zeros(1,i-1+(j-2)*m) -r zeros(1,m-1) t -1 zeros(1,m-2) -rzeros(1,(n-j)*m-i)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+left(y(j+1));end;%左下角的那个点i=1;j=1;a(1,:)=[t -1 zeros(1,m-2) -r zeros(1,(n-1)*m-1)];b(1)=h1*f(x(2),y(2))+r*bottom(x(2))+left(y(2));%右下角的那个点i=m;j=1;a(i+(j-1)*m,:)=[zeros(1,m-2) -1 t zeros(1,m-1) -r zeros(1,(n-2)*m)]; b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*bottom(x(i+1))+right(y(j+1)); %左上角的那个点i=1;j=n;a(i+(j-1)*m,:)=[zeros(1,(n-2)*m) -r zeros(1,m-1) t -1 zeros(1,m-2)]; b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*top(x(i+1))+left(y(j+1));%右上角的那个点i=m;j=n;a(i+(j-1)*m,:)=[zeros(1,(n-1)*m-1) -r zeros(1,m-2) -1 t];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*top(x(i+1))+right(y(j+1));u=a\bab2，精确解：function g=ni(x1,x2,y1,y2,M,N)m=M-1;n=N-1;x=x1+(x2-x1)*(0:M)/M;y=y1+(y2-y1)*(0:N)/N;for i=1:mfor j=1:nu1(i+(j-1)*m)=f1(x(i+1),y(j+1))endend（1）辅助函数function g=f(x,y)g=0;function g=bottom(x)g=2*log(x);function g=right(y)g=log(4+y^2);function g=top(x)g=log(x^2+1);function g=left(y)g=log(1+y^2);function g=f1(x,y)g=log(x^2+y^2);运行fivepoints(1,2,0,1,4,4)u =数值解0.4847467147016780.8376456266975491.1390195099193150.5944295076643080.9158860659528741.1974022894530100.7539416986884711.0340668399966291.287784599003526a =4 -1 0 -1 0 0 0 0 0 -1 4 -1 0 -1 0 0 0 0 0 -1 4 0 0 -1 0 0 0 -1 0 0 4 -1 0 -1 0 0 0 -1 0 -1 4 -1 0 -1 0 0 0 -1 0 -1 4 0 0 -1 0 0 0 -1 0 0 4 -1 0 0 0 0 0 -1 0 -1 4 -1 0 0 0 0 0 -1 0 -1 4b =0.5069117244448540.8109302162163292.5210301235267010.2231435513142101.4469189829363251.3872704470929461.1786549963416462.919669266564466运行ni(1,2,0,1,4,4)u1 =精确解Columns 1 through 30.485507815781701 0.838329190404443 1.139434283188365 Columns 4 through 60.594707107746693 0.916290731874155 1.197703191312341 Columns 7 through 90.753771802376380 1.034073767530539 1.287854288306638 误差很小（2）辅助函数function g=f(x,y)g=-4;function g=bottom(x)g=x^2;function g=right(y)g=(y-1)^2;function g=top(x)g=(x-2)^2;function g=left(y)g=y^2;function g=f1(x,y)g=(x-y)^2;fivepoints(1,2,0,1,4,4)fivepoints(0,1,0,2,4,4)u =0.062500000000000-0.0000000000000000.0625000000000000.5625000000000000.2500000000000000.0625000000000001.5625000000000001.0000000000000000.562500000000000a =Columns 1 through 32.500000000000000 -1.000000000000000 0 -1.000000000000000 2.500000000000000-1.0000000000000000 -1.000000000000000 2.500000000000000 -0.250000000000000 0 00 -0.250000000000000 00 0 -0.2500000000000000 0 00 0 00 0 0Columns 4 through 6-0.250000000000000 0 00 -0.250000000000000 00 0 -0.2500000000000002.500000000000000 -1.000000000000000 0 -1.000000000000000 2.500000000000000-1.0000000000000000 -1.000000000000000 2.500000000000000 -0.250000000000000 0 00 -0.250000000000000 00 0 -0.250000000000000Columns 7 through 90 0 00 0 00 0 0-0.250000000000000 0 00 -0.250000000000000 00 0 -0.2500000000000002.500000000000000 -1.000000000000000 0 -1.000000000000000 2.500000000000000-1.0000000000000000 -1.000000000000000 2.500000000000000b =0.015625000000000-0.1875000000000000.1406250000000000.750000000000000-0.250000000000000-0.2500000000000002.7656250000000000.3125000000000000.390625000000000精确解ni(0,1,0,2,4,4)u1 =u1 =Columns 1 through 30.062500000000000 0 0.062500000000000 Columns 4 through 60.562500000000000 0.2500000000000000.062500000000000Columns 7 through 91.562500000000000 1.0000000000000000.562500000000000误差很小（3）辅助函数function g=f(x,y)g=cosd(x+y)+cosd(x-y);function g=bottom(x)g=cosd(x);function g=right(y)g=-cosd(y);function g=top(x)g=0;function g=left(y)g=cosd(y);function g=f1(x,y)g=cosd(x)*cosd(y);数值解Pi=3.1415926fivepoints(0,pi,0,pi/2,4,4)u =0.6578183624886530.000000024999241-0.6578183271343870.5049807980892560.000000019229497-0.5049807708946410.2736443626241530.000000010432161-0.273644347870850a =10 -1 0 -4 0 0 0 0 0 -1 10 -1 0 -4 0 0 0 0 0 -1 10 0 0 -4 0 0 0 -4 0 0 10 -1 0 -4 0 0 0 -4 0 -1 10 -1 0 -4 0 0 0 -4 0 -1 10 0 0 -4 0 0 0 -4 0 0 10 -1 0 0 0 0 0 -4 0 -1 10 -1 0 0 0 0 0 -4 0 -1 10b =4.5582604075302670.000000137720159-4.5582602127645491.323957*********0.000000023374742-1.3239570281549570.7165204234523470.000000012650320-0.716520405562093精确解ni(0,pi,0,pi/2,4,4)u1 =Columns 1 through 30.653281493003155 0.000000024755257-0.653281457993935Columns 4 through 60.500000013397448 0.000000018946853-0.499999986602551Columns 7 through 90.270598066826879 0.000000010253963-0.270598052325585误差很小注：（1）需要对数值解与精确解作比较，以及不同步长选取下的误差比较。
泊松方程是在数学中的静电学，机械工程学和理论物理学中常见的偏微分方程。
它以法国数学家，几何学家和物理学家Poisson的名字命名。
泊松首先获得没有重力源的泊松方程△Φ= 0（即拉普拉斯方程）；考虑重力场时，△Φ= f（f为重力场的质量分布）。
后来，它扩展到了电场，磁场和热场分布。
该方程通常用格林函数法求解，但也可以用分离变量法和特征线法求解。
泊松方程为△φ=f
在这里△代表的是拉普拉斯算符(也就是哈密顿算符▽的平方)，而f 和φ 可以是在流形上的实数或复数值的方程。
当流形属于欧几里得空间，而拉普拉斯算子通常表示为，
因此泊松方程通常写成
在三维直角坐标系，可以写成
如果没有f，这个方程就会变成拉普拉斯方程△φ=0.
泊松方程可以用格林函数来求解;如何利用格林函数来解泊松方程可以参考screened Poisson equation[1] 。
现在有很多种数值解。
像是松弛法，不断回圈的代数法，就是一个例子。
数学上，泊松方程属于椭圆型方程(不含时线性方程)。
折叠编辑本段静电场的泊松方程
泊松方程是描述静电势函数V与其源(电荷)之间的关系的微分方程。
▽^2V=-ρ/ε
其中，ρ为体电荷密度(ρ=▽·D，D为电位移矢量。
)，ε为介电常
数绝对值εr*εo。
简称椭圆型方程,一类重要的偏微分方程。
早在1900年D.希尔伯特提的著名的23个问题中,就有三个问题(第19、20、23问题)是关于椭圆型方程与变分法的。
八十多年来，椭圆型方程的研究获得了丰硕的成果。
椭圆型方程在流体力学、弹性力学、电磁学、几何学和变分法中都有应用。
拉普拉斯方程是椭圆型方程最典型的特例。
拉普拉斯方程许多定常的物理过程，如稳定的热传导过程、牛顿引力理论及电磁理论中的位势、弹性薄膜的平衡、不可压流体的定常运动等，提出形如(1)的方程，称之为拉普拉斯方程，以及泊松方程(2)式中ρ一般有密度的意义。
容易得到方程(1)和(2)的一些特解。
由于方程是线性的，因此可以由已知的一些特解叠加而得到新的解。
积分也是一种叠加。
通过积分型叠加,便可得到方程(1)的如下的重要解：(3)式中S为一曲面，μ为定义在S上的连续函数。
由(3)确定的函数u在S以外的地方满足方程(1)。
非齐次方程(2)有一个重要的特解，它就是以ρ为密度的体位势：(4)只要ρ在域Ω内有界且连续可微,由(4)确定的函数u在Ω内就满足方程(2)，而在Ω外则满足方程(1)。
在应用上，往往不是求一些特解，而是求满足某些附加条件的解。
例如，第一边值问题（狄利克雷问题）：；第二边值问题（诺伊曼问题）：。
这里Ω为(x,y,z)空间的一个有界域，φ为定义在边界嬠Ω上的已知连续函数，n为嬠Ω的单位外法向量。
这些边值问题的解的惟一性，由调和函数的一个极值性质很容易推出。
拉普拉斯方程的二次连续可微解，称为调和函数。
极值原理域Ω内的调和函数不可能在域内一点取极大值或极小值，除非这个调和函数恒等于常数。
若调和函数的最大值只在某一边界点p上达到，则（假设u在p点可微）。
这些边值问题的解的存在性，也不难证明。
由格林公式可以推得(5)从而有式中，称为拉普拉斯方程关于域Ω的格林函数。
由此引出解第一边值问题的如下方法：先求出G(ξ,η,ξ；x,y,z)，再将所给的边界值代入，得到(6)只要对边界面再加上一些限制,就可以证明由(6)确定的函数u是第一边值问题的解。
