MATLAB求解PDE问题

MATLAB求解PDE问题（1）——概述、例子(转)(2011-07-2016:48:45)MATLABPDEToolbox提供利用有限元方法求解偏微分方程的GUI以及相应的命令行函数。利用该工具箱可以求解椭圆型方程、抛物型方程、双曲型方程、特征值方程以及非线性方程。PDEToolbox的功能非常强大，网上有许多利用PDEToolbox解决各种物理问题的论文，还有专门介绍工具箱的参考书。网上的例子虽然很多，但是大部分是介绍PDE工具箱自带的一些例子，这些例子中解的区域，边界条件是PDE工具箱已经编写好的，直接调用就可以。对于该如何自己设定求解区域及边界条件，却很少有人涉及。网上搜索发现只有刘平在博客中详细介绍过求解区域的设定。下面以一个椭圆型方程的例子来详细说明求解的各个步骤，希望对大家能有所帮助。设要求如下形式的椭圆方程的解：按照PDE的要求，将方程化为标准形式求解后的图像如下，第一幅图是解的图像，第二幅是计算误差。从第二幅图可以看到，计算的最大误差是10-3方量级。通过这个例子我们可以基本掌握PDE求解偏微分方程的步骤和方法，后面我将详细介绍如何设置区域及边界条件。掌握了区域和边界条件的设定，就可以轻松求解遇到的偏微分方程了。图后是附带的matlab命令以及注释，并提供m文件附件下载，下载后解压即可。希望能对大家有所帮助。下面是编写的求解上述方程的matlab语句及说明：g='mygeom';b='mybound';定义区域，边界条件。mygeom是定义区域的子函数名，函数名可根据自己的需要取定，区域的确定规则由pdegeom函数说明，注意pdegeom函数只是说明如何定义区域，它并不直接确定区域；mybound是定义边界条件的子函数名，与区域类似，边界的确定规则由函数pdebound确定。后面我会详细介绍区域和边界的取法。[p,e,t]=initmesh(g);网格初始化，此处也可以写成[p,e,t]=initmesh('mygeom');这样可以省略上面的语句[p,e,t]=refinemesh(g,p,e,t);[p,e,t]=refinemesh(g,p,e,t);加密网格两次，需要加密几次重复几次即可，根据具体问题确定加密次数U=assempde(b,p,e,t,1,0,'2*(x+y)-4');调用assempde函数计算方程的数值解，assempde函数的详细用法可以参考MATH网站或者PDE的使用指南。常用的用法是[u,res]=assempde(b,p,e,t,c,a,f)，其中b为边界条件，此处也可以写为'mybound'，p,e,t，为网格参数，c,a,f，为方程的参数，后面也可以加猜测值以及各种属性。pdesurf(p,t,U)gridon;xlabel('x');ylabel('y');zlabel('u')colorbarview([6030])画出解的图形。注意，为了让结果更直观一些，使用view函数调整了视点位置。大家可以自行调整视角，满意即可。exact=p(1,:).^2+p(2,:).^2-p(1,:).*p(2,:).*(p(1,:)+p(2,:));exact=exact';figurepdesurf(p,t,U-exact)gridonxlabel('x');ylabel('y');zlabel('error')colorbarview([6030])由于方程有解析解，我们可以比较数值计算的误差。如果能求得解析解，我们也不会设计各种方法求数值解了，因此，这一步在大多数情况下是用不上的，这里只是为了比较计算结果，验证计算的精度。MATLAB求解PDE问题（2）——确定几何区域(2012-06-1416:20:38)前一篇介绍了如何利用Matlab求解椭圆型方程，下面介绍如何确定求解的几何区域。PDEToolbox中规定几何区域的m文件是pdegeom.m。但是pdegeom并不是一个可以调用的函数，它只是规定了应该何如定义区域，具体的区域则要根据研究的问题来决定。函数pdegeom释义如下：参数为0个时，即没有参数时，返回边界的段数；参数为1个时，即只有bs，返回输出区域边界的参变量范围矩阵d；参数为2个时，返回每段边界长度为s时的坐标。函数参数意义bs表示指定的边缘线段，如矩形边界为四段，三角开边界肯定为三段…。s为第bs段线段弧长的近似（估计）值，bs与s可以为向量，但是要一一对应，即bs为几个值，s也得为几个值。输出变量[x,y]是每条线段起点和终点所对应的坐标。这个函数编制的关键是，函数内边界上的坐标((x(t),y(t))是用参变量t表示的，返回值是求得边界任意长度时的坐标（x(t),y(t)）值，参量可以有很多种选法。回到前一篇中，给定的方程的求解区域是[0,1,0,1]的一个正方形，我们将它命名为mygeom。下面我们来看下mygeom是怎么编写的。function[x,y]=mygeom(bs,s)nbs=4;%表示边界的段数ifnargin==0,x=nbs;%不给定输入变量时，输出表示几何区域边界的线段数returnendd=[0000%表示每条线段起始值的参量值t（四条边界，所以有四列）1111%表示每条线段的终点值参量值t0000%沿线段方向左边区域的标识值，如果是1，表示选定左边区域，如果%是0，表示不选定左边区域1111%沿线段方向右边区域的标识值，规则同上。];bs1=bs(:)';iffind(bs11|bs1nbs),error('PDE:squareg:InvalidBs','Nonexistentboundarysegmentnumber.')endifnargin==1,x=d(:,bs1);%给定一个输入变量时，输出区域边界数据的矩阵returnendx=zeros(size(s));y=zeros(size(s));[m,n]=size(bs);ifm==1&&n==1,bs=bs*ones(size(s));%expandbselseifm~=size(s,1)||n~=size(s,2),error('PDE:squareg:SizeBs','bsmustbescalarorofsamesizeass.');endif~isempty(s),%第一段边界ii=find(bs==1);iflength(ii)x(ii)=interp1([d(1,1),d(2,1)],[01],s(ii));%通过参量来确定边界上的值y(ii)=interp1([d(1,1),d(2,1)],[11],s(ii));%end%第二段边界ii=find(bs==2);iflength(ii)x(ii)=interp1([d(1,2),d(2,2)],[11],s(ii));y(ii)=interp1([d(1,2),d(2,2)],[10],s(ii));end%第三段边界ii=find(bs==3);iflength(ii)x(ii)=interp1([d(1,3),d(2,3)],[10],s(ii));y(ii)=interp1([d(1,3),d(2,3)],[00],s(ii));end%第四段边界ii=find(bs==4);iflength(ii)x(ii)=interp1([d(1,4),d(2,4)],[00],s(ii));y(ii)=interp1([d(1,4),d(2,4)],[01],s(ii));endend编辑完该函数后，可以利用pdeplot函数来测试设置的几何区域是否满足要求。在matlab命令窗口输入：pdegplot('mygeom'),axisequal[p,e,t]=initmesh('mygeom');pdemesh(p,e,t),axisequal结果如下：mygeom的m文件放在附件里，大家可以下载。本篇写作时参考了刘平的博客，已在前面指出，另参考了《偏微分方程的MATLAB解法》及《MATLABPDEToolboxUserGuide》，有兴趣的读者可以参阅。MATLAB求解PDE问题（3）——确定边界条件前一篇给出了如何确定几何区域，本篇接着给出如何确定边界条件，在几何区域和边界条件都确定好之后，就可以利用PDEToolbox对给定的PDE进行计算了。先回忆下前面的边界条件：上面的四个边界条件中，前两个是Dirichlet边界，后两个是Neumann边界。我们先了解下PDEToolbox中规定有三种边界条件，一是Dirichlet边界条件，一是广义Neumann条件，一是混合边界条件（只用于方程组）。这和传统的定义有所不同，按照传统的观点，PDEToolbox中的广义Neumann条件应该是混合边界条件（Robin边界条件），而传统的Neumann边界条件是指边界处的法向导数为零。我们先不考虑方程组，Dirichlet边界条件和广义Neumann条件写成如下形式：其中是外法线方向，是边界法向量与x轴的夹角。规定边界条件的m文件是pdebound函数，该函数并不提供边界条件，而是要求用户按照规定的方法给出边界条件，用户可以自己给定函数的名字。我们计算中给出的边界条件函数是mybound。调用格式为[q,g,h,r]=mybound(p,e,u,time)。pdebound函数中最主要的内容是bl矩阵，bl包括了所有的边界信息。bl的每一列都是边界条件矩阵的对应列，每一列都必须满足下面的规则：第一行是方程的维数N，一个方程N=1，两个方程构成的方程组，N=2，…；第二行是Dirichlet边界条件数M，Dirichlet边界条件时M=N，广义Neumann边界条件时M=0，如果是方程组，0MN之间表示混合边界条件；第三行到第3+N2-1行是表示q的字符串的长度，这个长度按与q有关的列方向的次序储存；第3+N2行到第3+N2+N-1是表示g的字符串的长度；第3+N2+N行到第3+N2+N+MN-1是表示h的字符串的长度，这个长度按与h有关的列方向的次序储存；第3+N2+N+MN行到第3+N2+N+MN+M-1是表示r的字符串的长度；接下来的行是包括MATLAB文本表达式所表示的真实边界条件函数。文本表达式是将实际的表达式通过double函数转化为ASCII码后的表达式。文本字符串的长度与前面四行规定的长度是一致的，两个字符串之间没有分隔符。可以插入含有下列变量的文本表达式：二维坐标x和y；边界线段参数s，弧长的比例。在边界线段的起始处s=0，向着线段的终点处增加到s=1；外法向量分量nx，ny。如果要用到切向量，可以用tx和ty表示，其中tx=-ny，ty=nx；解u（除非已指定输入参数u）；时间t（除非已指定输入参数time）。注意，在只有一个方程的时候，即N=1，如果M=0时，即是广义Neumann边界条件，此时不需要表示表示Dirichlet边界条件的h，r两行，因此表示q和g的字符串从第5行开始（前四行分别为N,M,q,g的指示行）；如果M=1，此时表示h，r的字符串从第9行开始（前六行分别是N,M,q,g,h,r的指示行，第七八行表示q，g，它们均为0，ASCII码为48）下面看下我们计算你的例子中边界mybound函数的内容：function[q,g,h,r]=mybound(p,e,u,time)%upper=double('2-2*x-x.^2');%y=1上边界条件函数%down=double('x.^2');%y=0下边界条件函数%left=double('y.^2');%x=0左边界条件函数%right=double('2-2*y-y.^2');%x=1右边界条件函数%upper=upper';%转化为列向量%down=down';%left=left';%right=right';%bl=[1111%N表示方程维数的指示行0011%M表示Dirichlet边界条件的指示行，M=0为广义Neumann边界条件1111%q表示q的字符串的长度，根据实际长度确定，10