Matlab解微分方程

1、第十六章 偏微分方程的数值解法科学研究和工程技术中的许多问题可建立偏微分方程的数学模型。包含多个自变量的微分方程称为偏微分方程(partial differential equation)，简称PDE。偏微分方程问题，其求解是十分困难的。除少数特殊情况外，绝大多数情况均难以求出精确解。因此，近似解法就显得更为重要。本章仅介绍求解各类典型偏微分方程定解问题的差分方法。16.1 几类偏微分方程的定解问题一个偏微分方程的表示通常如下： (16.1.1)式中，是常数，称为拟线性(quasilinear)数。通常，存在3种拟线性方程： 双曲型(hyperbolic)方程：； 抛物线型(parabolic

2、)方程：； 椭圆型(ellliptic)方程：。16.1.2 双曲型方程最简单形式为一阶双曲型方程： (16.1.2)物理中常见的一维振动与波动问题可用二阶波动方程： (16.1.3)描述，它是双曲型方程的典型形式。方程的初值问题为： (16.1.4)边界条件一般有三类，最简单的初边值问题为： (16.1.5)16.1.3 抛物型方程其最简单的形式为一维热传导方程： (16.1.8)方程可以有两种不同类型的定解问题：(1) 初值问题： (16.1.6)(2) 初边值问题： (16.1.7)其中，为已知函数，且满足连接条件： (16.1.8)边界条件为第一类边界条件。第二类和第三类边界条件为：

3、(16.1.9)其中，。当时，为第二类边界条件，否则称为第三类边界条件。16.1.4 椭圆型方程其最典型、最简单的形式是泊松(Poisson)方程 (16.1.10)特别地，当时，即为拉普拉斯(Laplace)方程，又称为调和方程： (16.1.11)Poisson方程的第一边值问题为： (16.1.12)其中为以为边界的有界区域，为分段光滑曲线，称为定解区域，分别为，上的已知连续函数。第二类和第三类边界条件可统一表示为： (16.1.13)其中为边界的外法线方向。当时为第二类边界条件，时为第三类边界条件。16.2 差分方法的基本概念差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的

4、数值解中应用最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点(网格点)集代替；将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组(称为差分格式)。如果差分格式有解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为原问题的近似解(数值解)。因此，用差分方法求偏微分方程定解问题一般需要解决以下问题：(1) 选取网格；(2) 对微分方程及定解条件选择差分近似，列出差分格式；(3) 求解差分格式；(4) 讨论差分格式解对于微

5、分方程解的收敛性及误差估计。下面，用一个简单的例子来说明用差分方法求解偏微分方程问题的一般过程及差分方法的基本概念。设有一阶双曲型方程初值问题。 (16.2.1)选取网格：图16.2.1 差分示意图首先对定解区域作网格剖分，最简单常用一种网格是用两族分别平行于轴与轴的等距直线： , (16.2.2)将分成许多小矩形区域。这些直线称为网格线，其交点称为网格点，也称为节点，和分别称作方向和方向的步长。这种网格称为矩形网格。(1) 对微分方程及定解条件选择差分近似，列出差分格式。如果用向前差商表示一阶偏导数，即： (16.2.3)其中。方程： (16.2.4)在节点处可表示为： (16.2.5)其中

6、：。由于当足够小时，在式中略去，就得到一个与方程相近似的差分方程： (16.2.6)此处，可看作是问题的解在节点处的近似值。同初值条件： (16.2.7)结合，就得到求问题的数值解的差分格式。式： (16.2.8)称为差分方程的截断误差。如果一个差分方程的截断误差为，则称差分方程对是阶精度，对是阶精度的。显然，截断误差的阶数越大，差分方程对微分方程的逼近越好。若网格步长趋于0时，差分方程的截断误差也趋于0，则称差分方程与相应的微分方程是相容的。这是用差分方法求解偏微分方程问题的必要条件。如果当网格步长趋于0时，差分格式的解收敛到相应微分方程定解问题的解，则称这种差分格式是收敛的。16.3 双曲

7、型方程的差分解法16.3.1 一阶双曲型方程的差分格式考虑一阶双曲型方程的初值问题： (16.3.1)将平面剖分成矩形网格，取方向步长为方向步长为，网格线为： , 为简便，记： , 以不同的差商近似偏导数，可以得到方程的不同的差分近似 (16.3.2) (16.3.3) (16.3.4)截断误差分别为，与。结合离散化的初始条件，可以得到几种简单的差分格式： (16.3.5) (16.3.6) (16.3.7)其中：。如果已知第层节点上的值，按上面三种格式就可求出第层上的值。因此，这三种格式都是显式格式。如果对采用向后差商，采用向前差商，则方程可化成： (16.3.8)相应的差分格式为： (16

8、.3.9)此差分格式是一种隐式格式，必须通过解方程组才能由第层节点上的值求出第层节点上的值。例16.3.1 对初值问题： 其中：，用差分格式求其数值解，取。解：记，由初始条件： 按差分格式： 计算公式为：。计算结果略。如果用差分格式： 求解，计算公式为： 计算结果略。与准确解： 比较知，按前一个差分格式所求得的数值解不收敛到初值问题的解，而后一个差分格式的解收敛到原问题的解。16.3.2 波动方程的差分格式对二阶波动方程： (16.3.1)如果令，则方程可化成一阶线性双曲型方程组： (16.3.2)记，则方程组可表成矩阵形式： (16.3.3)矩阵有两个不同的特征值，故存在非奇异矩阵，使得：

9、(16.3.4)作变换，方程组可化为： (16.3.5)方程组由二个独立的一阶双曲型方程联立而成。因此本节主要讨论一阶双曲型方程的差分解法。 下面给出如下波动方程和边界条件的差分格式： (16.3.6) (16.3.7) 将矩形划分成个小矩形，长宽分别为：，形成一个网格，如图16.3.1所示。图16.3.1 网格图可通过差分方程的方法计算近似值： 在连续行内，网格点的真实值为。 求和的中心差分为： (16.3.8) (16.3.9)在每一行的网格间距是均匀的：，而且。同时，它在每一列也是均匀的，。接下来，将和去掉，用(16.3.8)和(16.3.9)中的近似，并按顺序代入(16.3.6)式，可

10、得到差分方程： (16.3.10)可用它来近似(16.3.6)式。为方便起见，可将代入上式，可得： (16.3.11)设行和的近似值已知，可用上式求网格的行： (16.3.12)式中，。根据上式左边的4个已知值可得到近似值，如图16.3.2所示。图16.3.2 波动方程的空格样板 用上式时，必须注意，如果计算的某个阶段带来的误差最终会越来越小，则方法是稳定的。为了保证上式的稳定性，必然使。还存在其他一些差分方程方法以，称为隐格式法，它们更难实现，但对无限制。16.3.3 差分方法求解波动方程的MATLAB程序 求解区间，以(16.3.7)为边界条件的波动方程的差分方法程序。*function

11、U = finedif(f,g,a,b,c,n,m) %Input - f=u(x,0) as a string 'f'% - g=ut(x,0) as a string 'g'% - a and b right endpoints of 0,a and 0,b% - c the constant in the wave equation% - n and m number of grid points over 0,a and 0,b%Output - U solution matrix; analogous to Table 10.1%Initialize

12、parameters and U h = a/(n-1);k = b/(m-1);r = c*k/h;r2=r2;r22=r2/2;s1 = 1 - r2;s2 = 2 - 2*r2;U = zeros(n,m); %Comput first and second rows for i=2:n-1 U(i,1)=feval(f,h*(i-1); U(i,2)=s1*feval(f,h*(i-1)+k*feval(g,h*(i-1) . +r22*(feval(f,h*i)+feval(f,h*(i-2);end %Compute remaining rows of U for j=3:m, f

13、or i=2:(n-1), U(i,j) = s2*U(i,j-1)+r2*(U(i-1,j-1)+U(i+1,j-1)-U(i,j-2); endend U=U'*16.4 抛物型方程的差分解法以一维热传导方程： (16.4.1)为基本模型讨论适用于抛物型方程定解问题的几种差分格式。16.4.1 差分格式的建立首先对平面进行网格剖分。分别取为方向与方向的步长，用两族平行直线，将平面剖分成矩形网格，节点为：。为简便，记： ，。16.4.2 微分方程的差分近似在网格内点处，对分别采用向前、向后及中心差商公式： (16.4.2)一维热传导方程可分别表示为： (16.4.3) (16.4.4

14、) (16.4.5)由此得到一维热传导方程的不同差分近似： (16.4.6) (16.4.7) (16.4.8)上述差分方程所用到的节点各不相同。其截断误差分别为，和。因此，它们都与一维热传导方程相容。如果将式：中的用代替，则可得到又一种差分近似： (16.4.9)差分方程用到四个节点。由Taylor公式容易得出： (16.4.10)故其的截断误差为。因而不是对任意的，此差分方程都能逼近热传导方程：，仅当时，才成立。综上可知，用不同的差商公式可以得到微分方程的不同的差分近似。构造差分格式的关键在于使其具有相容性、收敛性和稳定性。前面三个方程都具有相容性，而此方程则要在一定条件下才具有相容性。1

15、6.4.3 初、边值条件的处理为用差分方法求解定解问题初值问题： (16.4.11)初边值问题： (16.4.12)还需对定解条件进行离散化。对初始条件及第一类边界条件，可直接得到： ， ， 式中：。对第二、三类边界条件： (16.4.13)需用差分近似。下面介绍两种较简单的处理方法。在左边界处用向前差商近似偏导数，在右边界处用向后差商近似，即： ， (16.4.14)则得边界条件的差分近似为： ， (16.4.15)其截断误差为。(2) 用中心差商近似,即： ， (16.4.16)则得边界条件的差分近似为： ， (16.4.17)其截断误差为。误差的阶数提高了，但出现定解区域外的节点和，这就

16、需要将解拓展到定解区域外。可以通过用内节点上的值插值求出和，也可以假定热传导方程在边界上也成立，将差分方程扩展到边界节点上，由此消去和。16.4.4 几种常用的差分格式以热传导方程的初边值问题： (16.4.18)为例给出几种常用的差分格式。(1) 古典显式格式令，则： (16.4.19)可改写成： 。将其与初始条件及第一类边界条件： ， ，结合，我们得到求解此问题的一种差分格式： (16.4.20)由于第0层上节点处的值已知，由此即可算出在第一层上节点处的近似值。重复使用此式，可以逐层计算出所有的，因此此差分格式称为典显式格式。又因式中只出现相邻两个时间层的节点，故此式是二层显式格式。(2)

17、 古典隐式格式将式： (16.4.21)整理并与初始条件及第一类边界条件式联立，得差分格式如下： (16.4.22)其中：。虽然第0层上的值仍为已知，但不能由上式直接计算以上各层节点上的值，必须通过解下列线性方程组： (16.4.23)式中，。 才能由计算，故此差分格式称为古典隐式格式。此方程组是三对角方程组，且系数矩阵严格对角占优，故解存在唯一。(3) Richardson格式Richardson格式是将式： 整理后与初始条件及第一类边界条件式联立。其计算公式为： 6.4.24)这种差分格式中所涉及的节点出现在三层上，故为三层显式格式。Richardson格式是一种完全不稳定的差分格式，因此

18、它在实际计算中是不能采用的。(4) 杜福特-弗兰克尔 (DoFort-Frankel) 格式DoFort-Frankel格式也是三层显式格式，它是由式： 与初始条件及第一类边界条件式结合得到的。具体形式如下： (16.4.25)用这种格式求解时，除了第0层上的值由初值条件得到，必须先用二层格式求出第1层上的值，然后再按上式逐层计算。(5) 六点隐式格式对二阶中心差商公式： 如果用在与处的二阶中心差商的平均值近似在处的值，即： 同时在点处的值也用中心差商： 近似，即： 这样又得到热传导方程的一种差分近似： (16.4.26)其截断误差为，将上式与初始条件及第一类边界条件式联立并整理，得差分格式：

19、 (14.5.27)此格式涉及到六个节点，它又是隐式格式，故称为六点隐式格式。与古典隐式格式类似，用六点格式由第层的值计算第层的值时，需求解三对角方程组： (16.4.28)式中，。此方程组的系数矩阵严格对角占优，故仍可用追赶法求解。例16.4.1 用古典显式格式求初边值问题。 的数值解，取。解：这里：，。由格式： 可得到： 将初值代入上式，即可算出：将边界条件，及上述结果代入又可求得：如此逐层计算，得全部节点上的数值解为：16.4.5 前向差分法求解热传导方程的MATLAB程序 设，其中，而且，其中，求解区间内的近似解。*function U=forwdif(f,c1,c2,a,b,c,n,

20、m) %Input - f=u(x,0) as a string 'f'% - c1=u(0,t) and c2=u(a,t)% - a and b right endpoints of 0,a and 0,b% - c the constant in the heat equation% - n and m number of grid points over 0,a and 0,b%Output - U solution matrix; analogous to Table 10.4%Initialize parameters and U h=a/(n-1);k=b/(m-

21、1);r=c2*k/h2;s=1-2*r;U=zeros(n,m); %Boundary conditions U(1,1:m)=c1;U(n,1:m)=c2; %Generate first rowU(2:n-1,1)=feval(f,h:h:(n-2)*h)' %Generate remaining rows of U for j=2:m for i=2:n-1 U(i,j)=s*U(i,j-1)+r*(U(i-1,j-1)+U(i+1,j-1); endend U=U'*16.4.6 Crank-Nicholson求解热传导方程的MATLAB程序 设，其中，而且，其中，求

22、解区间内的近似解。*function U=crnich(f,c1,c2,a,b,c,n,m) %Input - f=u(x,0) as a string 'f'% - c1=u(0,t) and c2=u(a,t)% - a and b right endpoints of 0,a and 0,b% - c the constant in the heat equation% - n and m number of grid points over 0,a and 0,b%Output - U solution matrix; analogous to Table 10.5%I

23、nitialize parameters and U h=a/(n-1);k=b/(m-1);r=c2*k/h2;s1=2+2/r;s2=2/r-2;U=zeros(n,m); %Boundary conditions U(1,1:m)=c1;U(n,1:m)=c2; %Generate first row U(2:n-1,1)=feval(f,h:h:(n-2)*h)' %Form the diagonal and off-diagonal elemnts of A and %the constant vector B and solve tridiagonal system AX=

24、B Vd(1,1:n)=s1*ones(1,n);Vd(1)=1;Vd(n)=1;Va=-ones(1,n-1);Va(n-1)=0;Vc=-ones(1,n-1);Vc(1)=0;Vb(1)=c1;Vb(n)=c2;for j=2:m for i=2:n-1 Vb(i)=U(i-1,j-1)+U(i+1,j-1)+s2*U(i,j-1); end X=trisys(Va,Vd,Vc,Vb); U(1:n,j)=X'end U=U'*16.5 椭圆型方程第一边值问题的差分解法本节以Poisson方程为基本模型讨论第一边值问题的差分方法。16.5.1 差分格式的建立考虑Poiss

25、on方程的第一边值问题： (16.5.1)图16.5.1取分别为方向和方向的步长，如图所示，以两族平行线：，将定解区域剖分成矩形网格。节点的全体记为： 。定解区域内部的节点称为内点，记内点集为。边界与网格线的交点称为边界点，边界点全体记为。与节点沿方向或方向只差一个步长的点和称为节点的相邻节点。如果一个内点的四个相邻节点均属于，称为正则内点，正内点的全体记为，至少有一个相邻节点不属于的内点称为非正则内点，非正则内点的全体记为。问题是要求出第一边值问题在全体内点上的数值解。 为简便，记： ，。对正则内点，由二阶中心差商公式： (16.5.2)Poisson方程在点处可表示为： (16.5.3)其

26、中： (16.5.4) 为其截断误差表示式，略去，即得与方程相近似的差分方程： (16.5.5)式中方程的个数等于正则内点的个数，而未知数则除了包含正则内点处解的近似值外，还包含一些非正则内点处的近似值，因而方程个数少于未知数个数。在非正则内点处Poisson方程的差分近似不能按上式给出，需要利用边界条件得到。边界条件的处理可以有各种方案，下面介绍较简单的两种。(1) 直接转移。用最接近非正则内点的边界点上的值作为该点上值的近似，这就是边界条件的直接转移。例如，点为非正则内点，其最接近的边界点为点，则有： (16.3.6)上式可以看作是用零次插值得到非正则内点处的近似值，容易求出，其截断误差为

27、。将上式代入，方程个数即与未知数个数相等。(2) 线性插值。这种方案是通过用同一条网格线上与点相邻的边界点与内点作线性插值得到非正则内点处值的近似。由点与的线性插值确定的近似值，得： (16.5.7)其中，其截断误差为。将其与方程相近似的差分方程联立，得到方程个数与未知数个数相等的方程组，求解此方程组可得Poisson方程第一边值问题的数值解。上面所给出的差分格式称为五点菱形格式， (16.5.8)实际计算时经常取，此时五点菱形格式可化为： (16.5.9)简记为： (16.5.10)其中：。例16.5.1 用五点菱形格式求解拉普拉斯(Laplace)方程第一边值问题。 其中：。取。解：网格中

28、有四个内点，均为正则内点。由五点菱形格式，得方程组： 代入边界条件： 其解为：，当时，对利用点，构造的差分格式，称为五点矩形格式，简记为： (16.5.11)其中：，其截断误差为： (16.5.12)五点菱形格式与矩形格式的截断误差均为，称它们具有二阶精度。如果用更多的点构造差分格式，其截断误差的阶数可以提高，如利用菱形格式及矩形格式所涉及的所有节点构造出的九点格式就是具有四阶精度的差分格式。16.5.2 用Dirichlet法求解Laplace方程的MATLAB程序 求解区间内的近似解。而且满足条件：，其中，而且，其中，而且存在整数和，使得：。*function U=dirich(f1,f2,f3,f4,a,b,h,tol,max1) %Input - f1,f2,