Ch7-求解常微分方程-1

1、225-232第七章求解偏微分方程常微分方程的求解是在高等数学里面讲过的，所以大家比较熟悉，对于偏微分方程的求解可能就有点陌生了。 事实上偏微分方程在工程上有着更广泛地应用，例如描述液体在多孔介质中的扩散，声学和电磁学中谐波的传播，热在固体中的传导等过程的微分方程都是偏微 分方程。所以求解偏微分方程在工程实际上有着非常重要的价值。手工求解常微分方程就已经非常麻烦和复杂了，那么求解偏微分方程就更加的复杂和繁MATLAB中有一个冗了。所幸的是有MATLAB这个计算工具来帮我们解决这一麻烦问题，专门用来求解片微分方程的工具箱PDE工具箱。这个工具箱不但提供有丰富的命令函数，使求解偏微分方程便的简单灵

2、活，而且该提供了一个求解偏微分方程的图形用户界面系统(GUI)，使得整个求解过程更加的人性化。特别是对于初学者来说 GUIjia ng更容易被接受,操作也更方便，所以偏微分方程的求解这一部分我们将主要介绍PDE工具箱的GUI系统。7-1偏微分方程的特点对于n阶常微分方程我们知道它的解取决于n个任意常数，例如一阶常微分方程的解可 以表示为：U = f(X)+ C ,c为任意常数。但是对于偏微分方程来说就不是由多少个常数来确 定的了，例如偏微分方程 - = f(x, y)的解可表示为:u(x,y)=隽 jy0f(jn)d上+w(x) +v(y)其中w(x)和v(y)为两个连续的任意可微函数。我们看

3、到偏微分方程的解可能是非常多 的，与常微分方程的解依赖于若干任意常数相比，它的自由度要大的多，对于多维偏微分方程的解更是这样。正是由于这个原因， 一般来说偏微分方程的解很难用通式表达出来。事实上我们常用到的是偏微分方程在某种特定条件下的解，这样靠着这些特定条件的约束我们就可以把偏微分方程的解表示出来了。我们把这些帮助确定偏微分方程特解的条件叫做定解条件，由于自变量在多维空间中的变化，其变化的区域非常复杂， 所以在区域边界上给出的定解条件也更加形式多样，我们一般称给定在区域边界上的定解条件为边界条件。我们看到边 界条件岁于求解 偏微分方程非常有用。由于偏微分方程的求解一般来说太过复杂，所以现在没

4、有一个对所有的偏微分进行求解的理论，所谓的求解偏微分方程也只是对于某些人们比较熟悉的类型的偏微分方程进行求rr 2j-2r 2CU2/C u , C u , C u、丰解。不同类型的偏微分方程代表不同的物理现象，例如方程=a (+ +)表ctex2cydZ示内部没有热源的情况下物体的温度分布应满足的条件,t表示时间变量,x、y、z分别表2 K示物体的三维空间坐标，这里 a2, K为物体的热传导系数，Q为该物体的热的容量。Q2 2 2同样t表再比如方程=a2（u+u）表示声波或者振动在介质中的传播，C texcycz即:椭圆 波动方程则 方程变为示振动传播的时间变量， X、y、z表示介质的三维空

5、间坐标，而 a则表示振动在该介质中传 播的速度。人们通过对偏微分方程的研究，把常见的偏微分方程归结为四种类型， 型方程、抛物型方程和混合型方程。前面讲过的热传导方程是一个抛物型方程， 是一个双曲型方程。如果热传导方程不考虑时间因素，很a2 （马+ =2 +雪）=0 ,这个方程就是一个椭圆方程了，混合型方程顾名思义就是一个泳dy辽方程具有上面三种类型方程中至少两种的特点的方程了。上面说过，边界条件对于求解偏微分方程很关键，在求解过程中最常用的地定解条件有两种，即：狄利克雷条件（又叫第一类边界条件）、诺曼条件（又叫第二类边界条件）。狄利克雷条件常写为udxj;对于波的振m = g（t, x ），即

6、满足给定区域边界Q的待求函数u为已知函数g;对于热传导方程狄利克雷条件的具体意义就是，已知物体的表面温度分布函数为g动方程则狄利克雷条件的具体意义就是， 波在介质的表面（或介质的给定区域的表面）动函数u为已知函数g。诺曼条件常写为 cn=®(t, xi)，Q(t, xi)为已知函数,Xi窃空是cn关于的外法线导数。对于热传导方程诺曼条件的具体意义为, 温度变化率为已知函数 0 法向的振动速率为已知函数7-2 PDE toolbox求解问题的背景知识沿给定区域Q边界的外法向的;对于波动方程则诺曼条件的具体意义为，沿给定区域Q边界的外0。1.求解方程的类型使用PDE工具箱求解偏微分方程的

7、基本类型有椭圆型方程、双曲型方程、抛物型方程、特征值方程、椭圆型方程组和非线形椭圆型方程组。PDE工具箱对这些方程的公式表述为：椭圆型方程:-V (cVu) +au = f,inO;双曲型方程:Q I Id -7 Ycu) +au = f ,in。抛物型方程:特征值方程:一可(c可U)+ au = kdu,inO这里可是微分算子，表示对多元函数求全微分运算。Q是给定的平面有界区域， 参数C、a、f和待求未知数u都是定义在Q上的函数。参数 d是定义在Q上的复函数，入是待求的特征值。在抛物型方程和双曲型方程中，参数函数C、a、f、d可以是时间t的函数。可以使用非线形解题器求解非线形椭圆型方程:-V

8、，（cNu） +a（u）u = f （u）,这里 c、a和f都是未知的待求函数 U的函数。另外，PDE工具箱提供的解题器都可以处理下面的方程组：T 可（5可5）可（彳可上）七115 +&灿2 = f1I一灯（56）可，（C22口2）+&215 +&22氏=f2利用命令还可以求解高阶方程组。2.边界条件PDE工具箱对它们的公式表我们知道常用的两个边界条件狄利克雷条件和诺曼条件, 示为：Dirichlet 条件：hu = r;Neumann 条件为： n "（c可u） +qu =g;这里n是/上的单位外法向矢量，h、r、q和g定义在60上的复值函数（对于特征 值问

9、题g=O,r=O ）。对于非线性问题，系数 g、q、h和r可以依赖于u;对于抛物型方程和双 曲型方程系数还可以依赖于时间 t。对于方程组情形，狄利克雷边界条件表示为：h11u1 +h12u2 =1215 中 h22u2 =r2般诺曼边界条件表示为:n（615）+ n （比帝上）+q11u1 +q12u2 =g1n ©2#5）+ n 亿22灯吐）+421山 +q22u2 = g?混合型边界条件表示为:h11u1 +02上=r1n '（C11 可uj + n '（02灯口2）+q11u1 +52口2 =g1 +01 n '（。21可山）+ n （c2u2）+q21

10、5 中q22u gh12这里计算卩时要注意满足狄利克雷条件。7-3 GUI求解偏微分方程的过程7-3-1 GUI求解问题的一般程序在命令窗口里输入 pdetool，回车，PDE工具箱的图形用户界面（GUI ）系统就启动了。整个过程大致可以分为六个阶从定义一个偏微分方程问题到完成解偏微分方程的定解，段，它们是：（1）绘制平面有界区域 Q ,被MATLAB称之为“ Draw模式”。系统提供了一系列的实体模型，包括矩形、圆、椭圆和多边形，我们可以通过公式把这些实体模型 组合起来，生成我们需要的平面区域定义边界，被MATLAB称之为“Boundary模式”。声明不同边界段的边界条件。定义偏微分方程，被

11、 MATLAB称之为“ PDE模式”。确定方程的类型和方程的 系数C、a、f和d。根据具体情况（例如材料的性质） ，我们还可以在不同的子 区域上独立的声明不同的系数。网格化Q区域，被MATLAB称之为“Mesh模式”。我们可以控制自动生成网格 的参数，还可以对生成的网格进行对次的细化，使网格分割的更细更合理。解偏微分方程，被 MATLAB称之为“ Solve模式”。对于椭圆型方程，我们可 以激活并控制非线性自适应解题器来处理非线性方程；对于抛物线方程和双曲线方程，设置初始边界条件后可以求出给定时刻t的解；对于特征值问题，可以求出给定区间上的特征值。求解完成后我们还可以返回到第四步，对网络进一步

12、细化，进行再次求解。（6）计算结果的可视化。 我我们可以通过设置系统提供的对话框，显示所求的解的表面图、网络图、等高线图和箭头梯度图（quiver）。对于抛物型和双曲型问题还可以进行动画演示。我们将在下一节中通过一个示如何通过上面的六个阶段来具体的求解一个偏微分方程, 例给出具体的操作步骤。7-3-2求解问题的示例让我们从求解最简单的偏微分方程一一泊松方程开始,泊松方程的表达式为:4川+卑+ %fex dy cz定解条件我们采用狄利克雷边界这里心是拉普拉斯算子，注意和微分算子V区分开。条件。我们看到泊松方程属于椭圆型偏微分方程，它可以表示内部没有热源的物体的热的分布，还可以表示静止的电磁场的分

13、布。首先在MATLAB的命令窗口输入 pdetool，回车确认，这样图形用户界面系统（ GUI） 就启动了。接下来我们将按照上一节所说的六个步骤开始操作了，启动后的GUI窗口如图7-1所示。R问E Tuulljuic - LDntitlcfdH回®FiLv EJh L atCh-dti± Urn-PDESdLvv PlLKaLp口 1 田.O|申 1 71 0Q PDE厶 血 1 =f jccnefc Scala X- 1J75Y: 3055Sei loihlUUc1D.6D.4D.20-0.2-D.il-0.608-1051 6Ini Cl-DiawiD Sjwniclr

14、fEKit图7-1 PDE工具箱的GUI窗口 在打开的窗口中，点击菜单Options，选中Grid选项，于是窗口中便出现了网格，这样方便我们在绘图时确定几何图形的位置和大小。用鼠标左键单击工具栏里的画圆按钮（左边第四个），使其处于选中状态（表现为下沉），然后在窗口的绘图区域按下左键拖动鼠标，于是一个椭圆就出现了，系统给它命名为E1。再选中工具栏里的画矩形的按钮（左边第二个），在绘图区域的中心位置拖动鼠标，于是一个矩形就出现了，系统给它命名为R1。把光标移到E1上面，双击鼠标左键，这时会弹出一个对话框， 用来定义E1的位置和大小，设置 如图7-2所示。丿 Object Dialog回区Objec

15、t type:EllipseX-center:-0.0549&1 S3206l0e88eY-center;Asemiaxes:073262442740091531B-semiawes:0.47175572519033975ElKRotation (degreec):Ndnrie：Cancel图7-2几何实体E1属性对话框然后选择0K按钮确认。用同样的办法调出R1的属性对话框，设置如下，然后确认,如图7-3所示。)Object Pialog叵区0 bject type：LeftBottom:Width:Name:RectangleOKCancel图7-3 R1属性对话框接下来把窗口中Se

16、t formula栏中的公式修改为 E1-R1，这样几何实体的绘制就完成了， 画完以后窗口如图 7-4所示。我们用鼠标点住画好的几何实体，然后拖动鼠标就可以移动实匸/jPPE ToolboxIDn±±1:lc<i|1E.4i t glL.niRnunduy PQE El«xhElifl tU<Lp ra壬很n PDEI=0* I Iciorwic ScaiffHlX： -aettnV: 1.9269StI larmiJac|Elfl1Irfo:A罰季t Tiw anckcfE prircird objecth). end pnrwmr abgt il

17、dwired.体的位置，另外还可以把画好的几何实体用 m文件的形式保存起来。图7-4绘好的几何实体EM曲，这样边界就生成了，如图接下来设置边界。用鼠标左键单击工具栏里的按钮7-5所示。Toplboj -in旦凶Vi ti JawC3| S|O| e| 3 gn| F口厶 血| = |7 |inwdc：匚皿Sol (ormUtacIl X: -1.003IMO Cbck and digd p«iirr«2 to ei«4轴日Iprt：.EhMf)Boiindai; Condition3回刃Bcundaiy condition equation:图7-5生成边界条件C

18、onditKT 痕:DescrphcnC Neunann0Co 訓 ci 已 nlValueDin畑orCancel我们看到边界是分段的并且有方向（带有箭头）。把光标移到一段边界上，然后双击鼠标左键，就会出现一个对话框，让我们设置该段的边界条件，如图7-6所示。本例中所有的边界都采用Dirichlet条件的默认值。图7-6边界条件对话框我们既可以采用上面的办法一段一段的设置边界条件，也可以一次设置好边界条件相同的边界段。先按下 shift键，然后用鼠标左键分别单击条件性同的边段，这样这些边界段都 被选中了，再使用上面的办法弹出对话框设置参数，于是被选中的段都被设置了。图7-8网格细化后的图像选

19、择方程的类型。单击工具栏里的PDE按钮，打开方程 类型对话框，由于 泊松方程适宜中椭圆型方程， 所以在对话框中选择 Elliptic 选项。方程的参数设置如图 7-7所示, 单击0K按钮确认。J PDE Specification匚叵冈'EgufliioriType ol PDE:CoelficientVdlUEf* Elipticc1协C Pflribolic4|c.o广 Hyperbo4cf|1D.O("Eigenmodes|协OKCancel图7-7方程类型对话框对设定的区域Q进行网格剖分。单击工具栏里的三角号按钮， 就可以进行网格剖分了。 我们还可以单击工具栏里网格剖分按钮右边的按钮（一个三角号里面还有一个小三角的按 钮）对网格剖分进一步细化，网格细化后的GUI窗口如图7-8所示。r-J F-DH riiL>l1jii>c 一 iP