一维偏微分方程的pdepe(matlab)函数解法

本文根据matlab帮助进行加工，根据matlab帮助上的例子，帮助更好的理解一维偏微分方程的pdepe函数解法，主要加工在于程序的注释上。ExamplesExample1.Thisexampleillustratesthestraightforwardformulation,computation,andplottingofthesolutionofasinglePDE.Thisequationholdsonanintervalfortimes.ThePDEsatisfiestheinitialconditionandboundaryconditionsItisconvenienttousesubfunctionstoplaceallthefunctionsrequiredbypdepeinasinglefunction.functionpdex1m=0;x=linspace(0,1,20);%linspace(x1,x2,N)linspace是Matlab中的一个指令，用于产生x1,x2之间的N点行矢量。%其中x1、x2、N分别为起始值、终止值、元素个数。若缺省N，默认点数为100t=linspace(0,2,5);sol=pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);%Extractthefirstsolutioncomponentasu.u=sol(:,:,1);%Asurfaceplotisoftenagoodwaytostudyasolution.surf(x,t,u)title('Numericalsolutioncomputedwith20meshpoints.')xlabel('Distancex')ylabel('Timet')%Asolutionprofilecanalsobeilluminating.figureplot(x,u(end,:))title('Solutionatt=2')xlabel('Distancex')ylabel('u(x,2)')%--------------------------------------------------------------function[c,f,s]=pdex1pde(x,t,u,DuDx)c=pi^2;f=DuDx;s=0;%--------------------------------------------------------------functionu0=pdex1ic(x)u0=sin(pi*x);%--------------------------------------------------------------function[pl,ql,pr,qr]=pdex1bc(xl,ul,xr,ur,t)pl=ul;ql=0;pr=pi*exp(-t);qr=1;Inthisexample,thePDE,initialcondition,andboundaryconditionsarecodedinsubfunctionspdex1pde,pdex1ic,andpdex1bc.Thesurfaceplotshowsthebehaviorofthesolution.Thefollowingplotshowsthesolutionprofileatthefinalvalueoft(i.e.,t=2).我们再将该问题复杂化，比如在原方程右边加一项2)(uuf，对于标准形式2us，其余条件不变functionpdex1m=0;x=linspace(0,1,20);%linspace(x1,x2,N)linspace是Matlab中的一个指令，用于产生x1,x2之间的N点行矢量。%其中x1、x2、N分别为起始值、终止值、元素个数。若缺省N，默认点数为100t=linspace(0,2,5);sol=pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);%Extractthefirstsolutioncomponentasu.u=sol(:,:,1);%Asurfaceplotisoftenagoodwaytostudyasolution.surf(x,t,u)title('Numericalsolutioncomputedwith20meshpoints.')xlabel('Distancex')ylabel('Timet')%Asolutionprofilecanalsobeilluminating.figureplot(x,u(end,:))title('Solutionatt=2')xlabel('Distancex')ylabel('u(x,2)')%--------------------------------------------------------------function[c,f,s]=pdex1pde(x,t,u,DuDx)c=pi^2;f=DuDx;s=u^2;%--------------------------------------------------------------functionu0=pdex1ic(x)u0=sin(pi*x);%--------------------------------------------------------------function[pl,ql,pr,qr]=pdex1bc(xl,ul,xr,ur,t)pl=ul;ql=0;pr=pi*exp(-t);qr=1;对比结果有差异，表示可行Example2.ThisexampleillustratesthesolutionofasystemofPDEs.Theproblemhasboundarylayersatbothendsoftheinterval.Thesolutionchangesrapidlyforsmall.ThePDEsarewhere.Thisequationholdsonanintervalfortimes.ThePDEsatisfiestheinitialconditionsandboundaryconditionsIntheformexpectedbypdepe,theequationsareTheboundaryconditionsonthepartialderivativesofhavetobewrittenintermsoftheflux.Intheformexpectedbypdepe,theleftboundaryconditionisandtherightboundaryconditionisThesolutionchangesrapidlyforsmall.Theprogramselectsthestepsizeintimetoresolvethissharpchange,buttoseethisbehaviorintheplots,theexamplemustselecttheoutputtimesaccordingly.Thereareboundarylayersinthesolutionatbothendsof[0,1],sotheexampleplacesmeshpointsnear0and1toresolvethesesharpchanges.Oftensomeexperimentationisneededtoselectameshthatrevealsthebehaviorofthesolution.程序在pdex4.m中functionpdex4%pdex4为文件名clcm=0;x=[00.0050.010.050.10.20.50.70.90.950.990.9951];t=[00.0050.010.050.10.511.52];sol=pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);%pdepe函数，用于直接求解偏微分方程，其形式为sol=pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)%下面为作图步骤u1=sol(:,:,1);%ui=sol(:,:,i)isanapproximationtotheithcomponentofthesolutionvectoru.u2=sol(:,:,2);figuresurf(x,t,u1)title('u1(x,t)')xlabel('Distancex')ylabel('Timet')figuresurf(x,t,u2)title('u2(x,t)')xlabel('Distancex')ylabel('Timet')%--------------------------------------------------------------function[c,f,s]=pdex4pde(x,t,u,DuDx)%被调用的函数，其作用是描述偏微分方程c=[1;1];f=[0.024;0.17].*DuDx;y=u(1)-u(2);F=exp(5.73*y)-exp(-11.47*y);s=[-F;F];%--------------------------------------------------------------functionu0=pdex4ic(x);%此函数是用来描述初值u0=[1;0];%--------------------------------------------------------------function[pl,ql,pr,qr]=pdex4bc(xl,ul,xr,ur,t)%此函数用来描述边界条件pl=[0;ul(2)];ql=[1;0];pr=[ur(1)-1;0];qr=[0;1];Inthisexample,thePDEs,initialconditions,andboundaryconditionsarecodedinsubfunctionspdex4pde,pdex4ic,andpdex4bc.Thesurfaceplotsshowthebehaviorofthesolutioncomponents.