《偏微分方程讲义建模、数值解和Matlab工具箱》由会员分享，可在线阅读，更多相关《偏微分方程讲义建模、数值解和Matlab工具箱（76页珍藏版）》请在金锄头文库上搜索。
1、偏微分方程讲义 建模、数值解和Matlab工具箱,提 纲,几个偏微分方程的典型问题建模 偏微分方程的基本概念 偏微分方程的建模 偏微分方程的数值方法 偏微分方程求解的数值方法 偏微分方程的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(
3、x)和T(x+ x)，利用力学基本定律，可以得到水平和竖直方向两个方程,这里，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到t
5、2流入该区域的热量为,因为,该区域的问题升高需要的热量满足,再根据能量守恒定理，传入的热量等于温度升高所需要的热量，于是有,后一方程进一步使用高斯公式,所以有,由t和V的任意性知道，被积函数恒为0，特别的，如果k是一个常量时，得到方程,当物体有内部热源的时候，方程为,因为,左边是升高温度需要的热量，右边第一项是传入的热量，第二项是内部热源的效果。,类似与弦振动方程，要确定温度，也需要初始条件和边界条件。,初始条件,第二边界条件，表示表面各点热流已知。,第一边界条件，表示表面各点温度已知。,第三边界条件，表示外界温度为u1,表面的热量和温度差成正比。,2.1 一些常见的偏微分方程,Poisson
6、 方程,带有稳定热源或内部无热源的稳定温度场的温度分布，不可压缩流体的稳定无旋流动及静电场的电势等均满足这类方程。下面的方程是Poisson 方程的第一边值问题。,第二第三边界条件,在研究热传导过程，气体扩散现象及电磁场的传播等随时间变化的非定常物理问题时，常常会遇到抛物型方程。下面是初边值问题,Cauchy问题,常见的一维振动与波动问题可用该方程表示。下面是初边值问题,双曲型方程,2.2 偏微分方程数值解,差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点（网格点）集代替；将问题
7、中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组（称为差分格式） 。如果差分格式有解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为原问题的近似解（数值解）。因此，用差分方法求偏微分方程定解问题一般需要解决以下问题： （i）选取网格； （ii）对微分方程及定解条件选择差分近似，列出差分格式； （iii）求解差分格式；,椭圆型方程第一边值问题的差分解法,首先将区域通过一系列横纵平行线划分，考察方程在交点处的值。注意到在边界的交点值是已知的，而在内部的交点的值是要求
8、的。为计算u在内部交点处的值，需要列出方程。使用二阶中心差商代替方程中的二阶导数，可以得到方程组,这里h和是横纵方向上的步长，右端是已知函数f在交点kj处的值，左边分子是待求的未知量。,一些内点的四个邻居都属于集合，这种点称为正则点，另一些内点的部分邻居不属于集合，这种点属于非正则点。对于每个正则点，可以得到一个方程，但是非正则点不能得到方程，这样最后得到的方程个数少于未知量个数。因此，对于非正则点要使用一些方法得到方程，常用的方法是直接转移和线性插值的方法。前者就是直接用内临近的点的值代替，后者是用内附近点的线性组合代替。 下面以例子讲解具体的求解过程。,例3，用五点菱形格式求解 Lapla
9、ce 方程第一边值问题,首先画出网格，如右图。可以看到，内点一共4个，都是正则点，其他的是边界点，其值已知，使用五点菱形格式，得到方程组如下：,代人边界条件的值，得到方程组,借助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)
10、=feval(f1,0:h:(n-1)*h); u(1:n,m)=feval(f2,0:h:(n-1)*h); b=-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
11、.2); f4=(y) log(4+y.2); a=1;b=1;h=1/3;tol=0.00001;max1=1000; sol=dirich(f1,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 %Outpu
12、t - U solution matrix; % If f1,f2,f3 and f4 are M-file functions call U = dirich(f1,f2,f3,f4,a,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
13、(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)=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 t
14、hroughout the grid -366- err=1; cnt=0;,while(errtol),上面的程序是使用SOR方法计算线性方程组。因为偏微分方程离散化之后得到的线性方程组比较特殊，同时规模非常大，所以往往选择比较有效的线性方程组数值算法。当问题规模较小，初学时可以不考虑这些。,抛物型方程的差分解法,使用差分进行近似，得到如下的方程,这里h和是横纵方向上的步长。,对初始条件及第一类边界条件，可直接得到,结合上页第一个差分格式，得到,计算的时候，由于第1层的数据已知，利用这些数据可以得到第二层的数据，再利用第二层的数据计算第三层，依次类推就能计算出全部的未知量。这个格式称为古典显
15、格式。,结合第二个差分格式，得到,这是一个三对角方程组，可以使用追赶法求解。这个格式称为古典隐格式。,必须通过解下列线性方程组,例4，用古典显示格式求初边值问题,的数值解，取h= 1, =0.5.,双曲型方程的差分解法,可以使用下面的差分进行近似，,上面的格式都是显格式。当然，也可以采用不同的格式组合得到其他的格式。这里不再讲解。,2.3偏微分方程的Matlab求解工具,PDEPE 命令的使用方法及例子,PDEPE（ initial-boundary value problems for parabolic-elliptic PDEs in 1-D ）是求解抛物线和椭圆型一维偏微分方程的命令。
16、为求解这些类型的各种各样的问题，在求解方程的形式，初始条件以及边界条件上，PDEPE使用下面的统一结构。方程形式,函数f用来描述“通量”，函数s用来描述“源”。,初始条件为,边界条件为（可以通过取不同的函数，对应全部的三类边界条件。）,Matlab的命令格式为pdepe。 M是参数，可以取0，1，2. pdefun是上面的函数f，icfun描述初始条件，bcfun描述边界条件，xmesh和tspan是x和t方向上的划分，options是选择项。下面的两个例子是matlab帮助文件中自带的例子，这里做一个详细的讲解。,例子5，求偏微分方程,其中x在0,1之间，且满足以下初边值条件,步骤1： 将欲求解的偏微分方程改写成如下的标准形式。,步骤 2: 编写偏微分方程的系数向量函数。 functio
