《《偏微分方程讲义》PPT课件.ppt》由会员分享，可在线阅读，更多相关《《偏微分方程讲义》PPT课件.ppt（76页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、偏微分方程讲义建模、数值解和Matlab工具箱温罗生2013.7提 纲几个偏微分方程的典型问题建模偏微分方程的基本概念偏微分方程的建模偏微分方程的数值方法偏微分方程求解的数值方法偏微分方程的Matlab求解工具Pdetool的使用方法练习题1.1 偏微分方程的基本概念 含有未知多元函数及其偏导数的方程，称之为偏微分方程。方程中出现的未知函数偏导数的最高阶数称为偏微分方程的阶。如果方程中对于未知函数和它的所有偏导数都是线性的，这样的方程称为线性偏微分方程，否则称它为非线性偏微分方程。一些常见的偏微分方程：拉普拉斯方程拉普拉斯方程：uxx+uyy+uzz=0适用于重力场或静电场。波波动动方程式方程
2、式：未知函数 u(x,y,z,t):热传导方程式：其中 k 代表该材料的热导率。泊松方程泊松方程：适用于所有物质或电荷的重力场或静电场。初始条件和边界条件称为定解条件，未附加定解条件的偏微分方程称为泛定方程。对于一个具体的问题，定解条件与泛定方程总是同时提出。定解条件与泛定方程作为一个整体，称为定解问题。1.2偏微分方程的建模例子1，弦的振动方程一根长为L的均匀细线，线密度为(x)，放于x轴上拉紧后，让它离开平衡位置做微小振动。求任意时刻t弦线的位移u(x,t)。使用元素法的思想，在弦上取x,x+x一段进行分析。因为弦的质量相对于张力非常小，因此可以看成弦仅仅受到两个张力，T(x)和T(x+x
3、)，利用力学基本定律，可以得到水平和竖直方向两个方程这里，a分别是两个力和水平方向的夹角，以及弦线在竖直方向的加速度。注意到弦仅仅在接近水平位置振动，所以和都是很小的量，于是前一个方程可以近似为这代表张力和位置无关。直观的理解这是显然的。同时因为和都是很小，也有如下的近似关系于是竖直方向的方程可以改写为注意到弧长ds近似等于dx，因此上述方程又可以写成当弦受竖直方向的外力F(x,t)作用时，竖直方向的方程可以写成因此方程可以写成仅仅通过上面的方程并不能完全确定弦的状态，我们必须先给出在初始时刻t=0时的状态，也就是所谓的初始时刻，即也就是给出初始位置和初始速度。根据方程和上面的初始条件，仍然不
4、能决定，还需要给出边界条件。更加具体的情况，一般有三类边界条件，第一类是两个端点固定，也就是第二类边界条件描述端点受到外力作用的情况，有第三类边界条件描述端点受到弹性支撑的情况，结合虎克定理，能得到条件例子2，热传导方程热传导定理的推导依据两个定理，一个是傅里叶定理和能量守恒定理。傅里叶定理说热量从高温区域向低温区域传播，热流强度满足所谓热流强度是单位时间通过单位横截面积的热量。上面dt是时间段，k是热传导系数。是温度场梯度方向面积的投影，u(x,y,z,t)是t时刻在空间点(x,y,z)的温度。取点(x,y,z)的一个小邻域，从时间t1到t2流入该区域的热量为因为该区域的问题升高需要的热量满
5、足再根据能量守恒定理，传入的热量等于温度升高所需要的热量，于是有后一方程进一步使用高斯公式所以有由t和V的任意性知道，被积函数恒为0，特别的，如果k是一个常量时，得到方程当物体有内部热源的时候，方程为因为左边是升高温度需要的热量，右边第一项是传入的热量，第二项是内部热源的效果。类似与弦振动方程，要确定温度，也需要初始条件和边界条件。初始条件第二边界条件，表示表面各点热流已知。第一边界条件，表示表面各点温度已知。第三边界条件，表示外界温度为u1,表面的热量和温度差成正比。2.1 一些常见的偏微分方程Poisson 方程带有稳定热源或内部无热源的稳定温度场的温度分布，不可压缩流体的稳定无旋流动及静
6、电场的电势等均满足这类方程。下面的方程是Poisson 方程的第一边值问题。第二第三边界条件在研究热传导过程，气体扩散现象及电磁场的传播等随时间变化的非定常物理问题时，常常会遇到抛物型方程。下面是初边值问题Cauchy问题常见的一维振动与波动问题可用该方程表示。下面是初边值问题双曲型方程2.2 偏微分方程数值解 差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连
7、续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。如果差分格式有解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为原问题的近似解（数值解）。因此，用差分方法求偏微分方程定解问题一般需要解决以下问题：（i）选取网格；（ii）对微分方程及定解条件选择差分近似，列出差分格式；（iii）求解差分格式；椭圆型方程第一边值问题的差分解法椭圆型方程第一边值问题的差分解法首先将区域通过一系列横纵平行线划分，考察方程在交点处的值。注意到在边界的交点值是已知的，而在内部的交点的值是要求的。为计算u在内部交点处的值，需要列出方程。使用二阶中心差商代替方程中的二阶导
8、数，可以得到方程组这里h和是横纵方向上的步长，右端是已知函数f在交点kj处的值，左边分子是待求的未知量。一些内点的四个邻居都属于集合，这种点称为正则点，另一些内点的部分邻居不属于集合，这种点属于非正则点。对于每个正则点，可以得到一个方程，但是非正则点不能得到方程，这样最后得到的方程个数少于未知量个数。因此，对于非正则点要使用一些方法得到方程，常用的方法是直接转移和线性插值的方法。前者就是直接用内临近的点的值代替，后者是用内附近点的线性组合代替。下面以例子讲解具体的求解过程。例3，用五点菱形格式求解 Laplace 方程第一边值问题首先画出网格，如右图。可以看到，内点一共4个，都是正则点，其他的
9、是边界点，其值已知，使用五点菱形格式，得到方程组如下：代人边界条件的值，得到方程组借助Matlab，可以轻松得到问题的解。下面是程序。clc,clear f1=(x)2*log(1+x);f2=(x)log(1+x).2+1);f3=(y)log(1+y.2);f4=(y)log(4+y.2);u=zeros(4);m=4;n=4;h=1/3;u(1,1:m)=feval(f3,0:h:(m-1)*h);u(n,1:m)=feval(f4,0:h:(m-1)*h);u(1:n,1)=feval(f1,0:h:(n-1)*h);u(1:n,m)=feval(f2,0:h:(n-1)*h);b=-
10、u(2,1)+u(1,2);u(4,2)+u(3,1);u(2,4)+u(1,3);u(3,4)+u(4,3);a=-4 1 1 0;1-4 0 1;1 0-4 1;0 1 1-4;x=ab当然，这里网格非常稀疏，得到结果比较粗糙，为得到更精确的值，可以将网格加密，比如两个方向各取1000个点。试自己修改上面的程序实现之。function sol=main f1=(x)2*log(1+x);f2=(x)log(1+x).2+1);f3=(y)log(1+y.2);f4=(y)log(4+y.2);a=1;b=1;h=1/3;tol=0.00001;max1=1000;sol=dirich(f1
11、,f2,f3,f4,a,b,h,tol,max1)%*function U=dirich(f1,f2,f3,f4,a,b,h,tol,max1)%Input-f1,f2,f3,f4 are boundary functions input as strings%-a and b right endpoints of 0,a and 0,b%-h step size%-tol is the tolerance%Output-U solution matrix;%If f1,f2,f3 and f4 are M-file functions call U=dirich(f1,f2,f3,f4,a,
12、b,h,tol,max1).%if f1,f2,f3 and f4 are anonymous functions call U=dirich(f1,f2,f3,f4,a,b,h,tol,max1).%Initialize parameters and U n=fix(a/h)+1;m=fix(b/h)+1;ave=(a*(feval(f1,0)+feval(f2,0).+b*(feval(f3,0)+feval(f4,0)/(2*a+2*b);U=ave*ones(n,m);%Boundary conditions U(1,1:m)=feval(f3,0:h:(m-1)*h);U(n,1:m
13、)=feval(f4,0:h:(m-1)*h);U(1:n,1)=feval(f1,0:h:(n-1)*h);U(1:n,m)=feval(f2,0:h:(n-1)*h);%逐次超松弛迭代法（SOR）参数 w=4/(2+sqrt(4-(cos(pi/(n-1)+cos(pi/(m-1)2);%Refine approximations and sweep operator throughout the grid-366-err=1;cnt=0;while(errtol)&(cnt 0.01,p,e,t=refinemesh(g,p,e,t);u=assempde(b,p,e,t,c,a,f);
14、%求得数值解 exact=(1-p(1,:).2-p(2,:).2)/4;err=norm(u-exact,inf);error=error err;end%结果显示 subplot(2,2,1),pdemesh(p,e,t);subplot(2,2,2),pdesurf(p,t,u)subplot(2,2,3),pdesurf(p,t,u-exact)例例9 考虑最小表面问题,在圆盘边界上x=u 2。这是椭圆型方程，其中g=circleg;b=circleb2;c=1./sqrt(1+ux.2+uy.2);rtol=1e-3;p,e,t=initmesh(g);p,e,t=refinemes
15、h(g,p,e,t);u=pdenonlin(b,p,e,t,c,0,0,Tol,rtol);pdesurf(p,t,u)例例10 求解正方形区域(x,y):-1 x,y 1上的热传导方程这里是抛物型方程，其中 d=1,f=0,a=0,c=1。%(1)问题定义 g=squareg;%定义正方形区域 b=squareb1;%边界上为零条件 c=1;a=0;f=0;d=1;%（2）产生初始的三角形网格 p,e,t=initmesh(g);%(3)定义初始条件 u0=zeros(size(p,2),1);ix=find(sqrt(p(1,:).2+p(2,:).2)0.4);u0(ix)=1%(4)
16、在时间段为0到0.1的20个点上求解 nframe=20;tlist=linspace(0,0.1,nframe);u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d);%(5)动画图示结果 for j=1:nframe pdesurf(p,t,u1(:,j);mv(j)=getframe;end movie(mv,10)例11,求解正方形区域(x,y):-1 x,y 1上的波方程这里是双曲型方程，其中d=1,f=0,a=0,c=1。%(1)问题定义 g=squareg;%定义正方形区域 b=squareb3;%定义边界 c=1;a=0;f=0;d=1;%（2）产生初始
17、的三角形网格 p,e,t=initmesh(g);%(3)定义初始条件 x=p(1,:);y=p(2,:);u0=atan(cos(pi*x);ut0=3*sin(pi*x).*exp(cos(pi*y);%(4)在时间段为0到5的31个点上求解 n=31;tlist=linspace(0,5,n);uu=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d);%(5)动画图示结果 for j=1:n pdesurf(p,t,uu(:,j);mv(j)=getframe;end movie(mv,10)2.4 pdetool的使用方法图形界面解法步骤大致上为：（1）定
18、义PDE 问题，包括二维空间范围，边界条件以及PDE 系数等。（2）产生离散化点，并将原PDE 方程离散化。（3）利用有限元素法(finite element method；FEM)求解并显示答案。前五个按钮为 PDE 系统的边界范围绘制功能，由左至右的用法为：1：以对角绘制矩形或正方形。按住鼠标左键可绘制矩形，而正方形需以按住右键的方式绘制。2：从中心点至某一角边的方式绘制矩形或正方形。同样地，鼠标左键绘矩形，右键绘正方形。3：由周围界线的方式绘制椭圆或圆形区域。鼠标左键用以绘制椭圆，而右键用来绘制圆形图形。4：以中心点向外的方式绘制椭圆或圆。同样地，鼠标左键及右键，分别用以绘制椭圆及圆形的
19、区域。5：用以绘制多边型等不规则区域，欲关闭此功能需按鼠标右键。在这些绘制按钮之后的7个按钮功能依序如下：6：用以给定边界条件。在此功能选定后，使用者可在任一图形边界上按住鼠标左键双击，然后在对话框中输入边界条件。7：用以指定 PDE 问题及相关参数。8：产生图形区域内离散化的网点。9：用以进一步将离散化的网点再取密一点(refine mesh)。10：在指定 PDE 系统，边界条件及区域后，按此钮即开始解题。11：用以指定显示结果绘制方式。12：放大缩小功能，便于图形绘制及显示。图形图形界面解法的使用步骤界面解法的使用步骤要利用pdetool 接口求解之前，需先定义PDE 问题，其包含三大部
20、份：（1）利用绘图(draw)模式，定义欲解问题的空间范围(domain)。（2）利用boundary 模式，指定边界条件。（3）利用PDE 模式，指定PDE 系数，即输入c，a，f 和d 等PDE 模式中的系数。在定义 PDE 问题之后，可依以下两个步骤求解（1）在mesh 模式下，产生mesh 点，以便将原问题离散化。（2）在solve 模式下，求解。（3）最后，在Plot 模式下，显示答案。以Poisson s方程式 u=10的求解为例，详细说明 pdetool 的用法。此问题的几何图形及相关边界条件，将于求解过程中加以说明。步骤 1：在命令窗口中键入pdetool 以进入GUI(gra
21、phical user interface)界面。选取Options 中之Grid 功能，以显示网格线。步骤 2：利用Draw 功能，画出问题之几何图形。请注意：使用者可利用内定对象“多边型“，“矩形”，“正方形”，“圆形”，及“椭圆型”，予以组合。步骤 3：选取 PDE 功能项，以输入 PDE 方程的系数及类型。因问题为 u=f，故此为椭圆型的问题，且其标准形式为(cu)+au=f，比较得知，c=1，a=0 和f=10，所以对话框输入的情况如图3。步骤 4：选取Boundary 功能，以输入边界条件。假设弧形部分边界条件为Neumann形，且为u n=5，与标准式比较知，g=-5 且q=0。
22、直线部分其边界条件则在 Dirichlettype 使h=0，r=0。对话框输入情况见图4。步骤 5：选取Mesh 功能，产生网点。使用者亦可进一步利用将网点取得密一点(refine mesh)，见图5。步骤 6：选取solve 功能，解此PDE，见图6。注意：1.MATLAB 会以图形的方式展示结果，使用者亦可点选plot 下的“parameters”功能，选择适当的方式显示图形及数据。2.另外，若使用者欲将结果输出到命令窗口中，以供后续处理，可利用solve 功能项下的“export solution”指定变量名称来完成。3.如果求抛物型或双曲型方程的数值解，还需要通过“solve”菜单下
23、的“parameters”选项设置初值条件。4.在上面定义边界条件和初始条件时，可以使用一些内置变量。（1）在边界条件输入框中，可以使用如下变量：二维坐标 x 和y，边界线段长度参数s，外法向矢量的分量nx 和ny（如果需要边界的切线方向，可以通过tx=-ny 和ty=nx 表示），解u。（2）在初值条件的输入框中，也可以输入用户定义的MATLAB 可接受变量（p，e，t，x，y）的函数。例 13 使用PDETOOL 重新求例10 的数值解。1）定义PDE 问题，包括二维空间范围，边界条件以及PDE 系数等。我们这里就省略了。2）区域剖分以后，通过“Mesh”菜单下的“Export Mesh”
24、选项可以把p，e，t三个参数分别输出到Matlab 工作空间。3）然后编写函数fun1(x,y)如下：function f=fun1(x,y);f=zeros(length(x),1);ix=find(x.2+y.20.16);f(ix)=1;其中的变量x，y 是MATLAB 可接受的内置变量。设置“solve”菜单下的“parameters”选项如下：时间框中输入：linspace(0,0.1,20)；初值框中输入：fun1。4）设置“plot菜单下的“parameters”选项如下：选择Height(3-D plot)和Animation两项。5）用鼠标点一下工具栏上的“”按钮，就可以画出
25、数值解的3-D 图形。例12,(MCM90A)药物在脑中的分布研究脑功能失调的人员要测试新的药物的效果，例如治疗帕金森症往脑部注射多巴胺（dopamine）的效果。为此，为了精确估计药物影响的脑部区域，他们必须顾及注射后药物在脑内分布区域的大小和形状。研究数据包括50各圆柱体组织样本的每个样本的药物含量的测定值（见附表和附图），每个圆柱体样本的长为0.76毫米，直径为0.66毫米，这些互相平行的圆柱体样本的中心位于网格距为1毫米*0.76毫米*1毫米的格点上，所以圆柱体互相间的底面上接触，侧面互不接触（见附图所示）。注射是在最高计数的那个圆柱体的中心附近进行的。自然在圆柱体之间以及由圆柱体样本
26、覆盖的区域外也有药物。试估计受到药物影响的区域中药物的分布。一个单元表示一个闪烁微粒的计数，或多巴胺的4.753*10-13克分子量，例如，附表指出位于后排当中的那个圆柱体的含药量是28353个单元。数学模型的研制还需要一些关于脑的生理生化特征及数据的附加假设。在研制模型的过程中，我们要证实这些假设中的合理性。例如，我们假定：可以粗略的认为大脑是均质的，扩散和衰减决定了脑中多巴胺的输运。我们还忽略了对流过程。此外，还假定了一次性注射，注射位置在y方向后排垂直柱体的中点。取样所需的确切时间不知道。由于注射到取样之间有足够的时间，因而取样所需时间相比之下可以忽略。决定多巴胺在大脑中的行为的物质可能极为复杂。通过假定分子扩散及组成元素衰减是脑中多巴胺输运的主要手段可得到合理的预测能力。所以我们假定该过程由一个二阶偏微分方程决定，这种形式的方程可应用于包括质量输运、流体动力学、热交换等多种问题。数值求积方法因为数据是用圆柱取样上的计数单位表示的，又因为模型预测了多巴胺的浓度值，有必要取应变量C在每个圆柱体上的积分值使之能和给定数据进行直接比较。3 练习题4.2011A题优秀论文以及MCM1999C题的优秀论文。