Ch7-求解常微分方程-1.doc

225-232第七章 求解偏微分方程常微分方程的求解是在高等数学里面讲过的，所以大家比较熟悉，对于偏微分方程的求解可能就有点陌生了。事实上偏微分方程在工程上有着更广泛地应用，例如描述液体在多孔介质中的扩散，声学和电磁学中谐波的传播，热在固体中的传导等过程的微分方程都是偏微分方程。所以求解偏微分方程在工程实际上有着非常重要的价值。手工求解常微分方程就已经非常麻烦和复杂了，那么求解偏微分方程就更加的复杂和繁冗了。所幸的是有MATLAB这个计算工具来帮我们解决这一麻烦问题，MATLAB中有一个专门用来求解片微分方程的工具箱PDE工具箱。这个工具箱不但提供有丰富的命令函数，使求解偏微分方程便的简单灵活，而且该提供了一个求解偏微分方程的图形用户界面系统（GUI），使得整个求解过程更加的人性化。特别是对于初学者来说GUIjiang更容易被接受，操作也更方便，所以偏微分方程的求解这一部分我们将主要介绍PDE工具箱的GUI系统。7-1 偏微分方程的特点对于n阶常微分方程我们知道它的解取决于n个任意常数，例如一阶常微分方程的解可以表示为：,c为任意常数。但是对于偏微分方程来说就不是由多少个常数来确定的了，例如偏微分方程的解可表示为： 其中w(x)和v(y)为两个连续的任意可微函数。我们看到偏微分方程的解可能是非常多的，与常微分方程的解依赖于若干任意常数相比，它的自由度要大的多，对于多维偏微分方程的解更是这样。正是由于这个原因，一般来说偏微分方程的解很难用通式表达出来。事实上我们常用到的是偏微分方程在某种特定条件下的解，这样靠着这些特定条件的约束我们就可以把偏微分方程的解表示出来了。我们把这些帮助确定偏微分方程特解的条件叫做定解条件，由于自变量在多维空间中的变化，其变化的区域非常复杂，所以在区域边界上给出的定解条件也更加形式多样，我们一般称给定在区域边界上的定解条件为边界条件。我们看到边界条件岁于求解偏微分方程非常有用。由于偏微分方程的求解一般来说太过复杂，所以现在没有一个对所有的偏微分进行求解的理论，所谓的求解偏微分方程也只是对于某些人们比较熟悉的类型的偏微分方程进行求解。不同类型的偏微分方程代表不同的物理现象，例如方程表示内部没有热源的情况下物体的温度分布应满足的条件，t表示时间变量，x、y、z分别表示物体的三维空间坐标，这里，K为物体的热传导系数，Q为该物体的热的容量。再比如方程表示声波或者振动在介质中的传播，同样t表示振动传播的时间变量，x、y、z表示介质的三维空间坐标，而a则表示振动在该介质中传播的速度。人们通过对偏微分方程的研究，把常见的偏微分方程归结为四种类型，即：椭圆型方程、抛物型方程和混合型方程。前面讲过的热传导方程是一个抛物型方程，波动方程则是一个双曲型方程。如果热传导方程不考虑时间因素，则方程变为，这个方程就是一个椭圆方程了，混合型方程顾名思义就是一个方程具有上面三种类型方程中至少两种的特点的方程了。上面说过，边界条件对于求解偏微分方程很关键，在求解过程中最常用的地定解条件有两种，即：狄利克雷条件（又叫第一类边界条件）、诺曼条件（又叫第二类边界条件）。狄利克雷条件常写为，即满足给定区域边界的待求函数u为已知函数g；对于热传导方程狄利克雷条件的具体意义就是，已知物体的表面温度分布函数为g ;对于波动方程则狄利克雷条件的具体意义就是，波在介质的表面（或介质的给定区域的表面）的振动函数u为已知函数g。诺曼条件常写为，为已知函数，是关于的外法线导数。对于热传导方程诺曼条件的具体意义为，沿给定区域边界的外法向的温度变化率为已知函数；对于波动方程则诺曼条件的具体意义为，沿给定区域边界的外法向的振动速率为已知函数。7-2 PDE toolbox求解问题的背景知识1求解方程的类型使用PDE工具箱求解偏微分方程的基本类型有椭圆型方程、双曲型方程、抛物型方程、特征值方程、椭圆型方程组和非线形椭圆型方程组。PDE工具箱对这些方程的公式表述为：椭圆型方程：双曲型方程：抛物型方程：特征值方程：这里是微分算子，表示对多元函数求全微分运算。是给定的平面有界区域，参数c、a、f和待求未知数u都是定义在上的函数。参数d是定义在上的复函数，是待求的特征值。在抛物型方程和双曲型方程中，参数函数c、a、f、d可以是时间t的函数。可以使用非线形解题器求解非线形椭圆型方程：这里c、a和f都是未知的待求函数u的函数。另外，PDE工具箱提供的解题器都可以处理下面的方程组：利用命令还可以求解高阶方程组。2边界条件我们知道常用的两个边界条件狄利克雷条件和诺曼条件，PDE工具箱对它们的公式表示为：Dirichlet条件：Neumann条件为：这里是上的单位外法向矢量，h、r、q和g定义在上的复值函数（对于特征值问题g=0,r=0）。对于非线性问题，系数g、q、h和 r可以依赖于u；对于抛物型方程和双曲型方程系数还可以依赖于时间t。对于方程组情形，狄利克雷边界条件表示为：一般诺曼边界条件表示为：混合型边界条件表示为：这里计算时要注意满足狄利克雷条件。7-3 GUI求解偏微分方程的过程7-3-1 GUI求解问题的一般程序在命令窗口里输入pdetool，回车，PDE工具箱的图形用户界面（GUI）系统就启动了。从定义一个偏微分方程问题到完成解偏微分方程的定解，整个过程大致可以分为六个阶段，它们是：（1） 绘制平面有界区域，被MATLAB称之为“Draw模式”。系统提供了一系列的实体模型，包括矩形、圆、椭圆和多边形，我们可以通过公式把这些实体模型组合起来，生成我们需要的平面区域。（2） 定义边界，被MATLAB称之为“Boundary模式”。声明不同边界段的边界条件。（3） 定义偏微分方程，被MATLAB称之为“PDE模式”。确定方程的类型和方程的系数c、a、f和d。根据具体情况（例如材料的性质）,我们还可以在不同的子区域上独立的声明不同的系数。（4） 网格化区域，被MATLAB称之为“Mesh模式”。我们可以控制自动生成网格的参数，还可以对生成的网格进行对次的细化，使网格分割的更细更合理。（5） 解偏微分方程，被MATLAB称之为“Solve模式”。对于椭圆型方程，我们可以激活并控制非线性自适应解题器来处理非线性方程；对于抛物线方程和双曲线方程，设置初始边界条件后可以求出给定时刻t的解；对于特征值问题，可以求出给定区间上的特征值。求解完成后我们还可以返回到第四步，对网络进一步细化，进行再次求解。（6）计算结果的可视化。我我们可以通过设置系统提供的对话框，显示所求的解的表面图、网络图、等高线图和箭头梯度图（quiver）。对于抛物型和双曲型问题还可以进行动画演示。如何通过上面的六个阶段来具体的求解一个偏微分方程，我们将在下一节中通过一个示例给出具体的操作步骤。7-3-2 求解问题的示例让我们从求解最简单的偏微分方程泊松方程开始，泊松方程的表达式为：这里是拉普拉斯算子，注意和微分算子区分开。定解条件我们采用狄利克雷边界条件。我们看到泊松方程属于椭圆型偏微分方程，它可以表示内部没有热源的物体的热的分布，还可以表示静止的电磁场的分布。首先在MATLAB的命令窗口输入pdetool，回车确认，这样图形用户界面系统（GUI）就启动了。接下来我们将按照上一节所说的六个步骤开始操作了，启动后的GUI窗口如图7-1所示。图7-1 PDE工具箱的GUI 窗口在打开的窗口中，点击菜单Options，选中Grid选项，于是窗口中便出现了网格，这样方便我们在绘图时确定几何图形的位置和大小。用鼠标左键单击工具栏里的画圆按钮（左边第四个），使其处于选中状态（表现为下沉），然后在窗口的绘图区域按下左键拖动鼠标，于是一个椭圆就出现了，系统给它命名为E1。再选中工具栏里的画矩形的按钮（左边第二个），在绘图区域的中心位置拖动鼠标，于是一个矩形就出现了，系统给它命名为R1。把光标移到E1上面，双击鼠标左键，这时会弹出一个对话框，用来定义E1的位置和大小，设置如图7-2所示。图7-2 几何实体E1属性对话框然后选择OK按钮确认。用同样的办法调出R1的属性对话框，设置如下，然后确认，如图7-3所示。图7-3 R1属性对话框接下来把窗口中Set formula栏中的公式修改为E1-R1，这样几何实体的绘制就完成了，画完以后窗口如图7-4所示。我们用鼠标点住画好的几何实体，然后拖动鼠标就可以移动实体的位置，另外还可以把画好的几何实体用m文件的形式保存起来。图7-4 绘好的几何实体接下来设置边界。用鼠标左键单击工具栏里的按钮，这样边界就生成了，如图7-5所示。图7-5 生成边界条件我们看到边界是分段的并且有方向（带有箭头）。把光标移到一段边界上，然后双击鼠标左键，就会出现一个对话框，让我们设置该段的边界条件，如图7-6所示。本例中所有的边界都采用Dirichlet条件的默认值。图7-6 边界条件对话框我们既可以采用上面的办法一段一段的设置边界条件，也可以一次设置好边界条件相同的边界段。先按下shift键，然后用鼠标左键分别单击条件性同的边段，这样这些边界段都被选中了，再使用上面的办法弹出对话框设置参数，于是被选中的段都被设置了。选择方程的类型。单击工具栏里的PDE 按钮，打开方程 类型对话框，由于 泊松方程适宜中椭圆型方程，所以在对话框中选择Elliptic选项。方程的参数设置如图7-7所示，单击OK按钮确认。图7-7 方程类型对话框对设定的区域进行网格剖分。单击工具栏里的三角号按钮，就可以进行网格剖分了。我们还可以单击工具栏里网格剖分按钮右边的按钮（一个三角号里面还有一个小三角的按钮）对网格剖分进一步细化，网格细化后的GUI窗口如图7-8所示。图7-8 网格细化后的图像求解方程。单击工具栏里的等号按钮，开始求解。这时窗口里就显示出了泊松方程在给定边界条件下的图形解，如图7-9所示。图7-9 泊松方程的图形解该图中的图像是通过颜色的深