泊松方程-椭圆型方程的五点格式
泊松方程的有限差分法的MATLAB实现
泊松方程的有限差分法的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 将上面的式子一次积分（注意符号）带入边界条件就能得出电场的分布，再次积分就能得出电势的分布。
泊松方程详解
椭圆型方程的五点差分格式
解 网格点如图所示
u(1,4)=180 u(2,4)=180 u(3,4)=180 u(0,3)=80
u(0,2)=80 U7 U4 U1 U8 U5 U2 U9 U6 U3 u(4,3)=0
u(4,2)=0
u(4,1)=0
u(0,1)=80
边界 x = 1上(3.15)
4U M ,m 2U M 1,m U M ,m1 U M ,m1 2hgM ,m
边界 y = 0上(3.16)
4U l ,0 2U l ,1 U l 1,0 U l 1,0 2hgl ,0
边界 y = 1上(3.17)
4U l , M 2U l , M 1 U l 1,M U l 1, M 2hgl , M
椭圆型方程的五点差分格式
例1 用五点差分格式求解 Laplace方程
u u 2 0 2 x y
2 2
在区域
( x, y) | 0 x 4,0 y 4
内的近似解，边界值为：
u( x,0) 20, u( x,4) 180,0 x 4
u(0, y) 80, u(4, y) 0,0 y 4
椭圆型方程的五点差分格式
§3.1 正方形区域中的Laplace 方程
Dirichlet边值问题的差分模拟
设Ω 是 xy 平面中的具有正方形边界 的 一个有界区域，考虑Laplace方程的第一边值 Dirichlet ）问题
2u 2u 2 2 0 ( x, y ) x y u ( x, y ) f ( x, y ) ( x, y )
U7 U8 U9 顶点
椭圆方程的公式
椭圆方程的公式椭圆方程是数学中一个非常重要的概念，它在物理、工程、计算机科学等领域都有广泛的应用。
本文将介绍椭圆方程的公式及其应用。
一、椭圆方程的定义椭圆方程是一个二元二次方程，其一般形式为：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足够小时，它近似为一个圆，因此可以用椭圆方程来描述圆形轨道。
数学物理方程泊松方程
在保险精算中，泊松方程可以用来预测未来的风险和 损失。
股票市场预测
在股票市场中，泊松方程可以用来预测股票价格的波 动和趋势。
泊松方程的扩展
04
非线性泊松方程
ห้องสมุดไป่ตู้
非线性泊松方程
在泊松方程的基础上，引入非线 性项，使其能够描述更复杂的物 理现象。
求解方法
由于非线性项的存在，求解非线 性泊松方程的难度增加，需要采 用迭代法、有限元法等数值解法。
泊松方程的来源和重要性
泊松方程的起源可以追溯到18世纪的数学和物理学领域。它是由法国数学家和物理学家西莫恩·德尼· 泊松在研究电场和重力场问题时提出的。
泊松方程在数学物理、工程技术和科学计算等领域具有广泛的应用价值。它涉及到许多物理现象和工 程问题的建模与求解，如静电场、位势论、量子力学和流体动力学等。因此，掌握泊松方程的基本理 论和方法对于深入理解和解决实际问题至关重要。
应用领域
非线性泊松方程在物理学、工程 学等领域有广泛的应用，如描述 晶体生长、流体动力学等。
泊松方程的数值解法
有限差分法
将泊松方程转化为差分方程，通过迭代求解。
有限元法
将求解区域划分为若干个小的单元，对每个单元进行近似求解，再 通过求解全局方程得到最终结果。
应用领域
数值解法广泛应用于实际问题的求解，如工程设计、物理模拟等。
泊松方程的应用
03
在物理中的应用
描述粒子在势场中的运动
泊松方程可以描述粒子在势场中的运 动，例如在量子力学和经典力学中， 粒子在势能场中的运动可以用泊松方 程来描述。
电磁波传播
热传导问题
在热传导问题中，泊松方程可以用来 描述温度场的变化和分布。
8_6泊松方程
Laplace方程的一般解
w( ρ , ϕ ) = C0 + D0 ln ρ + ∑ ρ m ( Am cos mϕ + Bm sin mϕ )
m =1
∞
+ ∑ ρ − m (Cm cos mϕ + Dm sin mϕ ).
m =1
8
∞
由于是圆内问题，因此在圆心处，解应为有限值。但是，在 一般解中lnρ 和 半径ρ负幂项在圆心处为无限大，所以应当排 除，即
10
小结 分离变量法、傅里叶级数法、冲量定理法和非其次边界条 件的处理方法。可以求解最一般的有界定解问题： 泛定方程－－非齐次
边界条件－－非齐次 初始条件－－非零值
11
一、一般的有界波动和输运问题 法1、边界条件齐次化。 泛定方程－－非齐次 泛定方程－－非齐次
边界条件－－非齐次 初始条件－－非零值
边界条件－－齐次 初始条件－－非零值
法2、利用叠加原理转化为两个简单的定解问题。
12
二、一般的有界稳定场问题 法1、用特解法，将非齐次方程转化为齐次方程。 法2、利用叠加原理转化为两个简单的可直接求解的定解问题。 注意: (1) 本章研究的全是定义在有界区域的定解问题，并且可 以用分离变数(傅立叶级数)法求解，但是并非任何有界的线性 的定解问题都能用此方法求解。例如变系数的线性偏微分方程
6
分离变量，令
w( ρ , ϕ ) = R( ρ )Φ(ϕ ),
Φ′′ + λΦ = 0,
得到两个常微分方程
ρ 2 R′′ + ρ R′ − λ R = 0,
利用周期性边界条件， Φ (ϕ ) = Φ (ϕ + 2π ), 得到本征值问题
五点差分格式
《微分方程数值解》大作业（一）——椭圆型方程编程计算：采用五点差分格式求如下椭圆型方程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）需要对数值解与精确解作比较，以及不同步长选取下的误差比较。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
x x( 2) x( 2 ) x( 2) d du ( p )dx (1) qudx (1) fdx, x x dx dx (1) ( 2) x( 2 ) x( 2)
或 W ( x ) W ( x ) (1) qudx (1) fdx, (2.14)
, 为给定常数.
区间剖分
将区间[a, b]分成N等分，分点为 xi a ih i 0,1,2, N , h (b a) / N . 于是我们得到区间I [a, b]的一个网格剖分， xi 称为网格结点（节点），间距h称为步长.
微分方程离散(差分方程）
现在将方程(1.1)在节点xi离散化，对充分光滑 的解u，由Taylor展式可得 u ( xi 1 ) 2u ( xi ) u ( xi 1 ) h2 d 2u ( x ) h 2 d 4u ( x ) [ ]i [ ]i O (h 3 ), (1.3) dx2 12 dx4 其中[ ]i 表示括号内函数在xi点取值.
第2章 poission方程的五点差 分法
§1 差分逼近的基本概念
考虑二阶常微分方程的边值问题 d 2u Lu 2 qu f a x b, (1.1) dx u ( a ) , u (b) (1.2) 其中q, f为[a, b]上的连续函数, q 0;
于是得到区间I的一个网格剖分，记hi xi xi 1 , 称h max hi为最大网格步长。用I h 表示网格内点
i
x1 , x2 , xN 1的集合, I h 表示内点和界点x0 a, xN b 的集合.
取相邻节点xi 1 , xi的中点x
i
1 2
1 ( xi 1 xi )(i 1,2, , N )， 2 x N b,
当h足够小，Ri (u )是h的二阶无穷小量. 若舍去Ri (u )，则得逼近方程(1.1)的差分方程： ui 1 2ui ui 1 Lhui qi ui f i , (1.6) 2 h 式中qi q ( xi ), f i f ( xi ).记[ Lu ]i f ( xi )， 称Ri (u )为差分方程(1.6)的截断误差.
i 1 2
p
u ( xi ) u ( xi 1 ) ] hi
ri [u ( xi 1 ) u ( xi 1 )] qi u ( xi ) hi hi 1
f i Ri (u ), (2.7)
1 d2 du 其中Ri (u ) (hi 1 hi )( [ 2 ( p )]i 4 dx dx 1 d 3u 1 d 2u [ p 3 ]i [r 2 ]i ) o(h 2 ), 12 dx 2 dx 为差分算子Lh的截断误差，舍去Ri (u ), 便得逼近边值问题(2.1), (2.2)的差分方程.
(2.22 )
也可用梯形公式，此时 2 pi 1 pi ai , pi 1 pi 1 d i (qi 1 qi 1 ), 2 2 2 1 f i ( f 1 f 1 ), i 2 i 2 2
(2.23)
§3 矩形网的差分格式
p( x
i
1) 2
u ( xi ) u ( xi 1 ) hi hi du d 3u [ p ] 1 [ p 3 ] 1 o(h 3 ), dx i 2 24 dx i 2 hi du d 3u [ p ] 1 [ p 3 ]i o(h 3 ), (2.4) dx i 2 24 dx
考虑Poisson方程 u f ( x, y ), ( x, y ) G (3.1) G是平面上一有界区域，其边界为分段光滑曲线， 在边界上满足下列边值条件之一：
u ( x, y ) (第一边值问题) (3.1)1 u ( x, y ) (第二边值问题) (3.1) 2 n u k u ( x, y ) (第三边值问题) (3.1)3 n
u0 , u N ,
有限体积法
考虑守恒型微分方程： d du ( p ) q ( x)u f ( x)， (2.13) dx dx 如果把它看作是分布在一根杆上的稳定温度场方 Lu 程，则在[a, b]内任一小区间 x (1) , x ( 2 ) ]上的热量守 [ 恒律具有形式
截断误差 ， (1.7) Ri (u )是用差分算子Lh 代替微分算子L所引起的 截断误差.
差分方程(1.6)当i 1,2, N 1, 时成立， 加上边值条件就得到关于ui的线性代数方程组： ui 1 2ui ui 1 Lhui qi ui f i , i 1,2, N 1, (1.8) 2 h u0 , u N . (1.9) 它的解ui是u ( x)于x xi的近似.称(1.8), (1.9)为逼近(1.1), (1.2)的差分方程或差分格式.此格式称为中心差分格式.
两族直线的交点(ih1 , jh2 )称为网点或节点， 记为( xi , y j )或(i, j ). 说两个节点( xi , y j )和( xi , y j )是相邻的, 如果 xi xi yi yi 1或 i i j j 1. h1 h2
2 2
p( x
i
1) 2
u ( xi 1 ) u ( xi ) hi 1 du h 2 i 1 d 3u [p ] 1 [ p 3 ]i o(h 3 ), (2.5) dx i 2 24 dx
由(2.5)减( 2.4)，并除以
hi hi 1 ，得 2 u ( xi 1 ) u ( xi ) u ( xi ) u ( xi 1 ) 2 [ p( x 1 ) p( x 1 ) ] i i hi hi 1 hi 1 hi 2 2
令p
1 i 2
p( x
1 i 2
), ri r ( xi ), qi q( xi ), f i f ( xi ),
则由(2.3), (2.6)知, 边值问题的解u ( x)满足方程：
Lhu ( xi ) u ( xi 1 ) u ( xi ) 2 [p 1 hi hi 1 i 2 hi 1
hi 1 hi 2 du du d 3u ([ p ] 1 [ p ] 1 ) [ p 3 ]i o(h 2 ), hi hi 1 dx i 2 dx i 2 12 dx hi 1 hi d 2 hi 1 hi d du du d 3u [ ( p )]i [ 2 ( p )]i [ p 3 )]i o(h 2 ), (2.6) dx dx 4 dx dx 12 dx
其中f ( x, y ), ( x, y ), ( x, y ), ( x, y )及k ( x, y ) 0都是 连续函数。
3.1 五点差分格式
取定沿x轴和y轴方向的步长h1和h2 , h (h h ) .
2 1 1 2 2 2
作两族与坐标轴平行的直线： x ih1 , i 0,1, y jh2 , j 0,1,
xi
利用中矩形公式，得 W ai [ 又 1 hi
i
1 2
ai
ui ui 1 , hi
(2.16) (2.17 ) (2.18)
xi
xi 1
dx 1 ] . p( x)
x
i
1 2
x
qudx
i
1 2
hi hi 1 d i ui , 2
将(2.16), (2.18)代到(2.14)，即得守恒型差分方程 [ai 1 ui 1 ui u u 1 ai i i 1 ] (hi hi 1 )d i ui , hi 1 hi 2 1 (hi hi 1 )i , 2 (2.20) (2.21)
i
1 2
1 2
,x
i
1 2
],
i
i
x
i
1 2
qudx
x
x
fdx, (2.14 )
i
1 2
考虑到p ( x)允许有间断点，由(2.15)进一步差分 化是不合适的。但“热流量” ( x)恒连续, W
故将(2.15)改写成 du W ( x) , dx p ( x) 再沿[ xi 1 , xi ]积分，得 ui ui 1 W ( x) dx, xi 1 p ( x )
于是在xi可将方程(1.1)写成 u ( xi 1 ) 2u ( xi ) u ( xi 1 ) q ( xi )u ( xi ) f ( xi ) Ri (u ), (1.4) 2 h h 2 d 4u ( x ) 其中 Ri (u ) [ ]i O(h 3 ), (1.5) 12 dx4
x x
其中 W ( x) p( x)
du ， (2.15) dx
把微分方程写成积分守恒型后，最高阶微商由二阶降到一阶， 从而可减弱对函数光滑性的要求。
特别于(2.14 )取[ x (1) , x ( 2 ) ]对偶单元[ x 则 W (x
1 ) W (x 2 1) 2 x
i 1 2
i
注意 : 方程(1.8)的个数等于网格内点x1 , x2 , , x N 1的 个数,因此它是N 1阶方程组.
§2 一维差分格式
考虑两点边值问题： d du du Lu ( p ) r qu f a x b, (2.1) dx dx dx u (a ) , u (b) ( 2 .2 ) 其中p C 1[a, b], p ( x) pmin 0, r , q, f C[a, b],
, 为给定常数.
我们将介绍差分格式的两种方法： 直接差分化法、有限体积法