matlab-pde选读

基于Matlab的偏微分方程数值解求数值解方法差分方法有限元方法MATLAB的pedpe函数MATLAB的PDEtool工具箱偏微分方程分类椭圆偏微分方程双曲线偏微分方程抛物线偏微分方程4椭圆偏微分方程特例—拉普拉斯方程拉普拉斯方程是最简单的椭圆偏微分方程，以下以拉普拉斯方程为例，讲述椭圆偏微分方程的的数值解法。拉普拉斯方程形式如下：5椭圆偏微分方程边界条件椭圆偏微分方程边界条件有以下三种提法：其中第一种提法最为普遍，下面以第一种边界条件，拉普拉斯方程为例介绍椭圆偏微分方程常用的五点差分格式和工字型差分格式的解法。6五点差分格式五点差分格式最常用的格式，其形式如下：注意：这里的边界为矩形区域。7五点差分格式算法注意：要保证x方向和y方向上的网格步长相等才能使用上面的公式。8五点差分格式在MATLAB中实现functionu=peEllip5(nx,minx,maxx,ny,miny,maxy)%x方向的节点数：nx%求解区间x的左端：minx%求解区间x的右端：maxx%y方向的节点数：ny%求解区间y的左端：miny%求解区间y的右端：maxy%求解区间上的数值解：uformatlong;hx=(maxx-minx)/(nx-1);hy=(maxy-miny)/(ny-1);u0=zeros(nx,ny);forj=1:nyu0(j,1)=EllIni2Uxl(minx,miny+(j-1)*hy);u0(j,nx)=EllIni2Uxr(maxx,miny+(j-1)*hy);endforj=1:nxu0(1,j)=EllIni2Uyl(minx+(j-1)*hx,miny);u0(ny,j)=EllIni2Uyr(minx+(j-1)*hx,maxy);end%边界条件的离散9五点差分格式在MATLAB中实现A=-4*eye((nx-2)*(ny-2),(nx-2)*(ny-2));b=zeros((nx-2)*(ny-2),1);fori=1:(nx-2)*(ny-2);ifmod(i,nx-2)==1ifi==1A(1,2)=1;A(1,nx-1)=1;b(1)=-u0(1,2)-u0(2,1);elseifi==(ny-3)*(nx-2)+1A(i,i+1)=1;A(i,i-nx+2)=1;%注意边界节点的离散方式b(i)=-u0(ny-1,1)-u0(ny,2);elseA(i,i+1)=1;A(i,i-nx+2)=1;A(i,i+nx-2)=1;b(i)=-u0(floor(i/(nx-2))+2,1);endendelseifmod(i,nx-2)==0ifi==nx-210五点差分格式在MATLAB中实现A(i,i-1)=1;%注意边界节点的离散方式A(i,i+nx-2)=1;b(i)=-u0(1,nx-1)-u0(2,nx);elseifi==(ny-2)*(nx-2)A(i,i-1)=1;A(i,i-nx+2)=1;b(i)=-u0(ny-1,nx)-u0(ny,nx-1);elseA(i,i-1)=1;A(i,i-nx+2)=1;A(i,i+nx-2)=1;b(i)=-u0(floor(i/(nx-2))+1,nx);endendelseifi1&&inx-2A(i,i-1)=1;A(i,i+nx-2)=1;A(i,i+1)=1;b(i)=-u0(1,i+1);elseifi(ny-3)*(nx-2)&&i(ny-2)*(nx-2)A(i,i-1)=1;A(i,i-nx+2)=1;%其余靠近边界点的离散A(i,i+1)=1;b(i)=-u0(ny,mod(i,(nx-2))+1);11五点差分格式在MATLAB中实现elseA(i,i-1)=1;A(i,i+1)=1;%与边界无关的点离散A(i,i+nx-2)=1;A(i,i-nx+2)=1;endendendendendu=A\b;%求线性方程组的解uformatshort12五点差分格式算例求解在命令窗口输入程序：u=peEllip5(51,0,2,51,0,2);13五点差分格式算例求解结果如果网格更密的话，即采用更多的节点进行计算，会得到更光滑的曲面。14工字型差分格式15工字型差分格式注意：这里给出的边界仍是矩形边界；并且做网格剖分时要保证x方向和y方向上的网格步长相等。16工字型差分格式给出同样的算例求解以下拉普拉斯方程：计算结果比较：用工字型差分格式计算得到的结果与五点差分格式得到的结果差不多，但是前者角点上的算法不太好，这是格式自身的缺陷。17双曲线偏微分方程——一维对流方程对流方程是最简单的双曲线偏微分方程，以下以一维、二维对流方程为例，讲述对流方程的数值解法。18一维对流方程——迎风格式其中a0,a0格式不同是为了满足差分格式的稳定性，若第一个式子a0,则差分格式是绝对不稳定的。上述迎风格式是条件稳定的并且是一阶精度的差分格式。19functionu=peHypbYF(a,dt,n,minx,maxx,M)%方程中的常数：a%时间步长：dt%空间节点个数：n%求解区间的左端：minx%求解区间的右端：maxx%时间步的个数：M%求解区间上的数值解：uformatlong;h=(maxx-minx)/(n-1);ifa0forj=1:(n+M)u0(j)=IniU(minx+(j-M-1)*h);%向左延拓M个节点的函数值endelseforj=1:(n+M)u0(j)=IniU(minx+(j-1)*h);%向左延拓M个节点的函数值endendu1=u0;fork=1:Mifa0fori=(k+1):n+Mu1(i)=-dt*a*(u0(i)-u0(i-1))/h+u0(i);endelsefori=1:n+M-ku1(i)=-dt*a*(u0(i+1)-u0(i))/h+u0(i);end一维对流方程——迎风格式20一维对流方程——迎风格式算例endu0=u1;endifa0u=u1((M+1):M+n);elseu=u1(1:n);endformatlong;21一维对流方程——迎风格式算例然后在MATLAB窗口输入下列命令：u=peHypbYF(1,0.005,101,0,1,100);22一维对流方程——迎风格式算例结果t=0时，u,x分布图t=0.5时，u,x分布图-0.5-0.4-0.3-0.2-0.100.10.20.30.40.500.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.60.723一维对流方程——拉克斯-弗里德里希斯格式24一维对流方程——拉克斯-弗里德里希斯格式25拉克斯-弗里德里希斯格式算例26拉克斯-弗里德里希斯格式算例27拉克斯-弗里德里希斯格式算例28一维对流方程——拉克斯-温德洛夫格式29一维对流方程——拉克斯-温德洛夫格式30一维对流方程——拉克斯-温德洛夫格式算例31一维对流方程——拉克斯-温德洛夫格式算例32一维对流方程——拉克斯-温德洛夫格式算例33一维对流方程——比姆-沃明格式34一维对流方程——比姆-沃明格式35一维对流方程——比姆-沃明格式算例36一维对流方程——多步格式多步格式也有多种，这里只简单介绍其中几种格式。包括Richtmyer多步格式、拉克斯-温德洛夫多步格式、MacCormack多步格式。37一维对流方程——多步格式38一维对流方程——多步格式算例Richtmyer多步格式算出的结果并不理想，不但左边有波动，而且光滑性也不好。拉克斯-温德洛夫多步格式算出的结果比较不错，虽然左边有点小波动，但是初始函数的宽度和高度都保持的不错。MacCormack多步格式求得的结果和拉克斯-温德洛夫多步格式算出的结果差不多。39双曲线偏微分方程——二维对流方程40二维对流方程——拉克斯-弗里德里希斯格式41二维对流方程——拉克斯-弗里德里希斯格式42二维对流方程——拉克斯-弗里德里希斯格式算例u=peHypb2LF(1,1,0.005,101,0,1,101,0,1,100);43二维对流方程——拉克斯-弗里德里希斯格式算例00.20.40.60.8100.5100.20.40.60.81结果与初始值对比，可以看出，拉克斯-弗里德里希斯格式算出的结果非常好。44二维对流方程——近似分裂格式近似分裂格式也是一种不错的格式，其结果也非常接近理论值。45抛物线偏微分方程——扩散方程在实际应用中遇到的抛物线偏微分方程主要是扩散方程。扩散方程有很强的物理背景，例如不用物质之间的扩散过程、热传递过程、波传播等过程都可以用扩散过程来描述。下面以扩散方程为例介绍几种差分格式。46扩散方程——显式格式47扩散方程——显式格式48扩散方程——显式格式算例u=peParabExp(1,0.005,101,0,1,100)49扩散方程——显式格式算例我们知道显示格式虽然简单，但其精度很差，而且求得的解容易出现震荡。次算例的结果如下：数值结果震荡的非常厉害，说明显式格式在这种条件下不稳定。00.10.20.30.40.50.60.70.80.91-8-6-4-202468x1021250扩散方程——跳点格式相同的算例数值结果同样震荡的非常厉害，说明跳点格式在这种条件下也不稳定。51扩散方程——隐式格式52扩散方程——隐式格式53扩散方程——隐式格式54扩散方程——隐式格式算例用隐式格式求解下面扩散方程的初值问题：u=peParabImp(1,0.005,101,0,1,0,1,100)55扩散方程——隐式格式算例结果是一条比较稳定的直线。00.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.60.70.80.9156扩散方程——克拉克-尼科尔森格式57扩散方程——克拉克-尼科尔森格式58扩散方程——克拉克-尼科尔森格式59扩散方程——克拉克-尼科尔森格式算例60扩散方程——克拉克-尼科尔森格式算例61扩散方程——加权隐式格式通过算例我们可以知道，通过取不同参数，用加权隐式格式算得的结果差别不大，只是达到稳态的时间不一样而已。62差分方法小结以上我们介绍了差分方法在椭圆型、抛物型和双曲型偏微分方程中的应用，用差分格式求解偏微分方程的基本步骤是一样的，首先把连续的问题离散化，建立差分格式，然后根据差分格式对求解区域进行网格剖分，最后求解方程。下面将简单介绍有限元方法。另外变分法、边界元法、混合有限元法和多重网格法等也是偏微分方程数值求解方法，有兴趣的同学可以参考相关书籍。63有限元方法——介绍有限元方法是数值求解偏微分方程边值问题的一种方法，此方法首先于20世纪50年代初由工程师提出，并用于求解简单的结构问题。有限元方法是这一种系统的数值方法，并奠定其数学基础，是在60年代中期以冯康先生为代表的中国学者与西方学者独立并行完成的。有限元方法不同于差分方法，主要有以下三大特点：(1)从数学物理问题的变分原理出发，而不是从微分方程出发，因此事从问题的整体描述而不是从问题的局部描述出发。(2)对所考虑问题的区域(以二维情形为例)作三角形（或其他简单多边形）剖分，而不是仅作矩形剖分。(3)用剖分区域上的简单函数(如分片多项式)去逼近原问题的解，而不是只在剖分节点上的数值逼近。64有限元方法——一维边值问题算例用有限元方法求解如下一维边值问题：22sin,012(0)0,u(1)0duxxdxu解：一维边值问题线性有限元数值解法的MATLAB程序如下：%线性有限元方法%参数设置N=10;X=0:1/(N+1):1;b=zeros(N+1,1);A=zeros(N+1,N+1);fori=2:N+1F1=@(x)2*(N+1)*(x-X(i-1)).*sin(pi*x/2);%句柄函数R1=quad(F1,X(i-1),X(i));F2=@(x)2*(N+1)*(X(i+1)-x).*sin(pi*x/2);R2=