第20章偏微分方程的数值解

/第二十章偏微分方程的 u

2u

f(x,

特别地，当f(x,y)0时，即为 斯

ux2y

2u2u x f(x, (x, u(x,y)

(x, (x,y为边界的有界区域，为分段光滑曲线，U称为定解区域，f(xy),(xy分别为上的已知连续函数。uu (x, (x,y 其中n为边界的外法线方向。当0时为第二类边界条件，0时为第三类边界 tax2 (a 初值问题（也称为Cauchy问题）

ax2

t x

u(x,0)

x

0tT,0xu(x,0)u(0,t)g(t),u(l,t)g 0t

其中(x),g1t),g2(t)为已知函数，且满足连(0)g1 (l)g2问题（7）中的边界条件u(0,tg1(tu(ltg2(tu

g 0t

u

g

0t其中1(t0,2(t0。当1(t2(t0时，为第二类边界条件，否则称为第三uau

2ax

2 2ax t0,xu(x,0) xx

2ax

t0,0x

tt

0xu(0,t)g

u(l,t)g2 0t §2考虑Poisson方程的第一边值问题（3）2u2u x f(x, (x,u(x,y)

(x,y

(x, h,分别x方向y方向的步长，以两族平行xxkkhyyjj(k,j0,1,2,L)将定解区域剖分成矩形网格。节点的全体记为Rxkyj|xkkhyjji,j为整数。定解区域内部的节点称为内点，记内点集RI为h。边界与网格线的交点称为边界点，边界点全体记为h。与节点(xk,yj)沿x方向或y方向只差一个步长的点(xk1,yj)和(xk,yj1)称为节点(xk,yj)的相邻节点。如果一个内点的四个相邻节点均属于U，称为正则内点，正则内点的全体记为(1)，至少有一个相邻节点不属于U的内点称为非正则内点，非正则内点的全体记为(2) 为简便记，记(kj)xkyju(k,j)u(xkyj),fkjf(xkyj(kj(1)u(k1,j)2u(k,j)u(k1,j)(k,ju(k,j1)2u(k,j)u(k,j1)O2(k,j Poisson方程（1）在点(k,ju(k1,j)2u(k,j)u(k1,j)u(k,j1)2u(k,j)u(k,jfk,jO(h22

在式（12）中略去O(h22，即得与方程（1）uk1,j2uk,juk1,

uk,j12uk,juk,j1

fk,

式（13）中方程的个数等于正则内点的个数，而未知数uk,j则除了包含正则内点处解u的近似值，还包含一些非正则内点处un直接转线性插

1

uk1,

uk,j

uk,j

k,

)

k,

1h2uk,jfk, 其中uk,juk1,juk1,jukj1ukj14ukj例 用五点菱形格式求解Lace方程第一边值问2u2u x (x,u(x,y)

lg[(1x)2y2

(x,y 其中{(x,y)|0x,y1}。取h 3当h时，利用点(k,jk1,j1k1,j1)12h2(uk1,j1uk1,j1uk1,j1uk1,j14uk,j)fk, 11

uk,

fk,

ukjuk1,j1uk1,j1uk1,j1uk1,j14ukj a 0 (a 首先对xt平面进行网格剖分。分别取h,x方向与t方向的步长，用两族平行xxkkh(k0,1,2,Lttjjj0,1,2,Lxt平面剖分成矩形格，节点为(xktj)(k0,1,2,L,j0,1,2,L)。为简便起见，记(k,j)(xk,yju(k,j)u(xk,yj)，k(xk)，g1jg1(tj)，g2jg2(tj)，1j1(tj)2j2(tj在网格内点(k,j

，对x2采用u(k,j1)u(k,j)au(k1,j)2u(k,j)u(k1,j)O(h2 u(k,j)u(k,j1)au(k1,j)2u(k,j)u(k1,j)O(h2 hu(k,j1)u(k,j1)au(k1,j)2u(k,j)u(k1,j)O(h2) uk,j1uk,

auk1,j2uk,juk1,jh

uk,juk,j1auk1,j2uk,juk1,j huk,j1uk,j1auk1,j2uk,juk1,j uk,0u(xk,0) (k0,1,L或k0,1,L, u0,ju(0,tj)un,

u(l,

j)g2

(j0,1,L,m n

lmT 在左边界(x0)x，在右边界(xl (n,j

u(1,j)u(0,j)hu(n,j)u(n1,j)h

(j0,1,L,u1,ju0,j h 1j0, 1h

(j0,1,L, un,jun1,j

2jn, 2 (0,(n,j

u(1,j)u(1,j)O(h2u(n1,j)u(n1,j)O(h2

(j0,1,L,u1,ju1,j

1j0, 1

(j0,1,L, n1,

2

n,

g2)和(n1,u值插值求出u1,j和un1,j，也可以假定热传导方程（5）在边界上也成立，将差分方程扩展到边界节点上，由此消去u1,j和un1,j。(i)古典显式格为便于计算，令r uk,j1ruk1,j(12r)uk,jruk1,（21（22）uk,j1ruk1,j(12r)uk,jruk1, (k1,2,L,n1,j0,1,L,muuk

(k1,2,L, u0,jg1j un,jg2 (j1,2,L,0层j0)上节点处的u值已知(uk,0k，由式（25）即可算出uj1上节点处的近似值uk,1。重复使用式（25），可以逐层计算出各层节点的近似值（21（22）uk,0 (k0,1,L, u0,jg1j un,jg2 (j0,1,L,其中r 。虽然第0层上的u值仍为已知，但不能由式（30）直接计算以上各层点上的值uk,j，故差分格式（26）称为古典隐式格式杜福特—克尔（DoFort—Frankel）格 1uk,j112r(uk1,juk1,j)12ruk,j1(k1,2,L,n1,j1,2,L,m (k0,1,L,

g1j un,jg2 (j0,1,L,求出第1层上的值uk,1，然后再按格式（27）逐层计算ukj(j2,3,Lm)。如果令v1

2axtv2 x，则方程（10）

a2v

2 记v(v1,v2 ，则方程组（28）可表成矩阵形v

a2 0xA PAP1 00 0wPvw1w2)T，方程组（29）w uau t x x长为，网格xxkkh(k0,1,2,Lttjjj0,1,2,L。为简见，记(k,j)(xkyju(k,j)u(xkyjk(xkuk,juk,j1uk, uk1,juk,auk,j1uk,auk,juk1,uk,j1uk,auk1,juk1, 一维状态空间的偏微分方程的解c(x,t,u,u)uxm(xmf(x,t,u,u))s(x,t,u,u x 其中时间介于t0ttfx则介于[ab有限区域之间。m式中f(xtuux)一项为流通量(flux)，而s(xtuux)为来源(source)项。角线元素一定不全为0。偏微分方程初始值可表示为

u(x,t0)v0p(x,t,u)q(x,t)f(x,t,u,ux)

x为两端点位置，即a或 命令pdepe的用法solpdepe(m,pdepe,icfun,bcfun,xmesh,tspan,x0a(起点)xNb(终点)tspan为时间变量t的向量，即tspant0t1,LtM，其中t0为起始时间，tM为c,f,spdefun(x,t,u,亦即，使用者仅需提供偏微分方程中的系数向量。cfs均为行(column)向量，而向量c即为矩阵c的对角线元素。icfun提供解u的起始值，其格式为uicfun(x值得注意的是upl,ql,pr,qlbcfunxl,ul,xr,ur,t其中ul和ur分别表示左边界(xlb和右边界(xrauplqlpqprqrpq的行向量。sol为解答输出。sol为的输出向量，sol(:,:i)为ui的输出，即uisol(:,:,i)。元素ui(jk)sol(jki)表示在ttspan(j)xxmesh(k时ui之答案。optionsodeset的使用方法。PDE求解器pdepe的算法，主要是将原来的椭圆型和拋物线型偏微分mesh点，以二阶空间离散化(spatialdiscretization)技术为之(KeelandBerzins,1990)ode15s的指令ode15sode解法，主要是因为在离散化的过程中，椭圆型偏微分方程被以ode15s便可顺利求解。x的取点(mesh)位置对解的精确度影响很大，若pdepe求解器给出“…hasmesh点数。另外，若状态u在某些特定点上有较快速的变动时，亦需将此处的点取密集些，以增加精确度。值得注意的是pdepe并不会自动做xmesh的自动取3以上。tspan的选取主要是基于使用者对那些特定时间的状态有而选定。而间(stepsize)的控制由程序自动完成uout,duoutdxpdeval(m,xmesh,ui,m代表问题的对称性。m=0表示平板；m=1表示圆柱体；m=2pdepe中的自变量m时间固定为tjtspan(j)下的解。duoutdxdudxref.Keel,R.D.andM.Berzins,“AMethodfortheSpatialDiscritizationofParabolicEquationsinOneSpaceariable，SIAMJ.Sci.andSat.Comput.,Vol.11,pp.1-例 2u 其中0x1起始值条IC：u(x,0)sin(BC1：u(0,t)BC2：etu(1,t)注：本问题的解析解为u(xt)etex20_1.m，然后求解。步骤 将欲求解的偏微u程改写成如式的标准式 x0x0u x x cx,t,u,ux和m0

fx,t,u,uxsx,t,u,ux0步骤 function3编写起始值条件。步骤4 （37， BC1：u(0,t)0

(0,t) xBC2：et1u(1,t) 即

xplu(0,和pret

qlqr步骤5 x=linspace(0,1,20);%x取20点步骤6 利用pdepe求解。m=0;%依步骤1之结果步骤7 ylabel('时间间为2(终点)时，各位置下的解，可输入以下指令(利用pdeval指令)：figure(2);%绘成图2xout=linspace(0,1,100);%输出点位置plot(xout,uout);%绘图综合以上各步骤，可写成一个程序求解 2。其参考程序如下functionu=sol(:,:,1);%取出答案xlabel('位置x')ylabel('tylabel('tM=length(t);%取终点时间的下表xout=linspace(0,1,100);%输出点位置plot(xout,uout);%绘图例

F(uu2 u2

2u

F(uu2 F(u1u2exp(5.73(u1u2exp(11.46(u1u2，且0x1和t0。初值条边值条

u1(x,0)u2(x,0)u1(0,t)u2(0,t)u1(1,t)u2(1,t)解步骤1：改写偏微分方程 u11*u1 uxF(u1u2 tu

x

F

)

2c

2x

u1 f u 2 xsF(u1u2F(uu) 01 xu 2

2 0.170 即0 pl ，ql 2

x同理，右边界条件x1处)为

u11 u 1 0 u 即

x pru11，qr 步骤2：编写偏微分方程的系数向量函数。c=[11]';f=[0.024s=[-FF]';步骤3：编写初始条件函数u0=[10]';pl=[0ul(2)]';ql=[1pr=[ur(1)-1qr=[0x=[00.0050.010.050.990.995t=[00.0050.010.050.10.511.5functionx=[00.0050.010.050.990.995t=[00.0050.010.050.10.511.5u1=sol(:,:,1);%第一个状态之数值解输出%%pdefunction[c,f,s]=ex20_2pdefun(x,t,u,dudx)c=[1f=[0.024s=[-FF]';functionu0=ex20_2ic(x)u0=[1function[pl,ql,pr,qr]=ex20_2bc(xl,ul,xr,ur,t)ql=[10]';pr=[ur(1)-1qr=[04 1T r(H e A p GCr rr p D2 1f r e A ar0 ur rr 0L0,TT0(r),ff0r0,Tf rrw,kerhw(TTwr

,f反应速率3H 环己1kK3KP3rA(1KH

BHKPKP)B ClnK12100RT32.3lnKH15500RT31.9RlnKB11200RT23.1RlnKC8900RT19.4R操作条件及物性数总 Pt

rw2.5cmG631kgm30

m2

y0Mav PB1200kg CP1.74kcalkg Hr49250kcalkg

h065.8kcalT(0)125Lu8.03mhrKe0.65kcal

m2hr0mhr0

hw112D0.755e

m2hr0率f表示，方能求解偏微分方程。基于以下的反应方程3H 2P m3HP

t1m3fmf t1m3fPP t1m3程序设计T r T rAB(Hr) 0L 1 r f r 1 rr f A

3L

r T rAB(Hr)p r pc，f ，s r f A r 和m1(圆柱)。另外，左边界条件r0处)01* f 即 pl，ql 同理右边界条件rrw)hw(TTw GCp

*f 即

prhw(TTw)，qrGCp 根据以上的分析，可编写程序求解此PDE问题，其参考程序如下functionglobalPtrwTwGMy0Mavrho_BCpdHrh0uRkehwPt=1.25;%总压(atm)rw=0.025;%管径(m)Tw=100+273;%壁温(℃)G=631;%质量流率(kg/m2hr)T=sol(:,:,1f=sol(:,:,2%PDEfunction[c1,f1,s1]=ex20_3_1pdefun(r,L,u1,DuDr)globalPtrwTwGMy0Mavrho_BCpdHrh0uRkehwDe%%%%c1=[1f1=[ke/(G*Cp)functionu0=ex20_3_1ic(x)u0=[125+2730]';functionglobalPtrwTwGMy0Mavrho_BCpdHrh0uRkehwDepl=[00]';ql=[1qr=[G*Cp例 B液体接触面之浓度为CA00.01molm3而改变，即在操作时间内(h10天)维持定值。气体A在液体B

2109

ABA与BABC，rA其反应速率常数k2107s1

气体B 1Bc3A在液体BCA CA(z,0)0，z

AzA

(0,t)

，t0

0

t

A

z

N

与标准式（35）比较，可得C1fDABCAzs0m0左边界(z0)plCA(0,t)CA0ql0。右边界zL)pr0qr

与标准式（35）m0C1

DAB

zskCA利用以上的处理结果，可编写参考程序如下functionglobalDABk%casefori=1:length(t)title('case(a)')xlabel('lengthylabel('time(day)')zlabel('conc.(mol/m^3)')xlabel('time(day)')ylabel('flux(mol/m^2.day)')%casefori=1:length(t)%title('case(b)')xlabel('lengthylabel('time(day)')zlabel('conc.(mol/m^3)')xlabel('time(day)')ylabel('flux(mol/m^2.day)')PDE%caseglobalDABkCA0%casefunction[c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)globalDABkCA0functionCA_i=ex20_3_2ic(z)globalDABkCA0§4二维状态空间的偏微分方程的解法下面我们讨论的方程是定义在平面上的有界区域 上，区域的边界记作 工具箱可以解决下列类型的偏微分方程：(cu)au in其中ca,f和未知的u可以是du(cu)au其中ca,fd可以依赖于时间t

indt2(cu)au in(cu)au in其中d是(c(u)u)a(u)u其中ca,f可以是u

f in(c11u1)(c12u2)a11u1a12u2(c21u)(c

)a

au

21

22 （i）Dirichlet条件：hu ann条件：n(cu)qu rh11u1h12u2r1h21u1h22u2r2annrn(c21u1)n(c22u2)q21u1q22u2（iii）n(c21u1)n(c22u2)q21u1q22u2g2这里的计算是使得满足Dirichlet例 2u求解区域为单位圆盘，边界条件为在圆盘边界上u0。解它的精确解为u(x,y)1x2y4 error=[];err=1;whileerr>0.01,u=assempde(b,p,e,t,c,a,f);%求得数值解1111|u u) {(x,y)|x2y2在圆盘边界上ux2解这是椭圆型方程，其中c

,a0,

1|u例8求解正方形区域{(xy|1xy1}u初始条件为u(0)

x2y2边界条件为Dirichlet条件u0解这里是抛物型方程，其中c1,a0,g='squareg';%b='squareb1'

0,d1for例9求解正方形区域{(xy|1xy1}2u 初始条件为u(0)arctan(cos(x))， x1上满足Dirichlet条件u0y1

ann条件 0解这里是双曲型方程，其中c1,a0,f0,d1