《最新[解微分方程]解微分方程的步骤.doc》由会员分享，可在线阅读，更多相关《最新[解微分方程]解微分方程的步骤.doc（65页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、最新解微分方程解微分方程的步骤微分方程求解求解微分方程 ：简单地说，就是去微分（去掉导数），将方程化成自变量与因变量关系的方程（没有导数）。近来做毕业设计遇到微分方程问题，搞懂后，特发此文，来帮广大同学，网友。1.最简单的例子：1.1 dy=2x y=x2+C dxdy1.2 求微分方程 =2xy的通解。 dx解 方程是可分离变量的，分离变量后得dy=2xdx y两端积分 ： dy, y=2xdx得： lny=x2+C1,y=ex2从而 ：又因为 +C1=eC1ex2。 eC1仍是任意常数，可以记作C2y=Cex。1.3 非齐次线性方程2y 求方程y-=(x+1)2的通解. x+1解：非齐次线
2、性方程。先求对应的齐次方程的通解。 5dy2y-=0， dxx+1dy2dx=， yx+1y=C(x+1)2用常数变易法：把C换成u(x)，即令y=u(x+1)2 （1）则有 dy=u(x+1)2+2u(x+1)， dx12， 代入原方程式中得 u=(x+1)两端积分，得 2u=(x+1)2+C。 33再代入（1）式即得所求方程通解2y=(x+1)2(x+1)2+C。 3法二： 假设待求的微分方程是:我们可以直接应用下式-P(x)dxP(x)dxy=e(Q(x)edx+C) 3dy+P(x)y=Q(x) dx得到方程的通解，其中，2， Q(x)=(x+1)2 P(x)=-x+1代入积分同样可得
3、方程通解 52y=(x+1)(x+1)2+C， 323 2.微分方程的相关概念：(看完后你会懂得各类微分方程)一阶微分方程：y=f(x,y)或P(x,y)dx+Q(x,y)dy=0可分离变量的微分方程：一阶微分方程可以化为g(y)dy=f(x)dx的形式，解法：g(y)dy=f(x)dx得：G(y)=F(x)+C称为隐式通解。dyy=f(x,y)=j(x,y)，即写成的函数，解法：dxxydydududxduy设u=，则=u+x，u+=j(u)，=分离变量，积分后将代替u，xdxdxdxxj(u)-ux齐次方程：一阶微分方程可以写成即得齐次方程通解。一阶线性微分方程：dy1+P(x)y=Q(x
4、)dx-P(x)dx当Q(x)=0时,为齐次方程，y=CeP(x)dx-P(x)dx当Q(x)0时，为非齐次方程，y=(Q(x)edx+C)edy2+P(x)y=Q(x)yn，(n0,1)dx全微分方程：如果P(x,y)dx+Q(x,y)dy=0中左端是某函数的全微分方程，即：uudu(x,y)=P(x,y)dx+Q(x,y)dy=0=P(x,y)=Q(x,y) xyu(x,y)=C应该是该全微分方程的通解。二阶微分方程：f(x)0时为齐次d2ydy+P(x)+Q(x)y=f(x) 2dxdxf(x)0时为非齐次二阶常系数齐次线性微分方程及其解法：(*)y+py+qy=0，其中p,q为常数；求
5、解步骤：1、写出特征方程：(D)r2+pr+q=0，其中r2，r的系数及常数项恰好是(*)式中y,y,y的系数；2、求出(D)式的两个根r1,r23、根据r1,r2的不同情况，按下表写出(*)式的通解：y+py+qy=f(x)，p,q为常数f(x)=elxPm(x)型，l为常数；f(x)=elxPl(x)coswx+Pn(x)sinwx型3.工程中的解法：四阶定步长Runge-Kutta算法 其中 h 为计算步长，在实际应用中该步长是一个常数，这样由四阶 Runge-Kutta算法可以由当前状态变量Xt 的值求解出下状态变量Xt +1 的 值 亲们，你们满意吗？ 一阶微分方程的解 一阶微分方程
6、的常数变易法的应用探析 The exploration of linear ordinary differential equation of first order with method of leadingvariables 作 者：刘 *专 业：数学与应用数学指导老师：杜 * *完成时间：2010年9月1号 摘 要常数变易法是作为求解一阶线性方程的解法给出的。本文先介绍一阶线性非齐次微分方程的常数变易法，然后讨论四种形式的一阶非线性微分方程的常数变易法，包括齐次方程、贝努力方程和黎卡提方程等的常数变易法。Method of leading variables is method of
7、solving linear ordinary differential equation of first order. This paper first introduces first-order differential equations of nonhomogeneous linear method, and then discuss variation of four types of first order nonlinear differential equation of variation, including homogeneous equation, the baye
8、sian equation and li CARDS carry equation of variation law.关键词：一阶线性； 一阶非线性； 常数变易法Key words: A linear ; First-order nonlinear ; Method of leading variables 目 录 1、一阶线性非齐次微分方程的常数变易法 . 4 2、一阶非线性微分方程的常数变易法 . 5dyyy=+g()xx . 5 2.1 齐次方程dxdy=p(x)y+Q(x)yn2.2 贝努力方程：dx . 5 dy=p(x)y2+Q(x)y+R(x)2.3 黎卡提方程：dx . 6yy
9、+p(x)e=Q(x)的微分方程 . 8 2.4 形如 目前，由于常微分方程应用的广泛性，人们基本满足于各类型方程的各自求解方法。基于此，常微分方程课程可以说是各类型的孤立技巧与方法的汇编，从内容联系上势必感到松散。因此，把握解常微分方程的方法，在学习此类课程时，不仅仅是记住一些解法，更重要的是强调思维方法的训练。由于常数变易法是求解微分方程的一种很重要的方法,且常应用于一阶线性微分方程的求解，因此本文对这一部分的内容做一系统整理。在数变易法中，将常数换成u(x)就可以得到非齐次线性方程的通解。1、一阶线性非齐次微分方程的常数变易法dy=p(x)y+Q(x) （1） dxdy=p(x)y (2
10、) 先解对应的其次线性微分方程dx为求解一阶非齐次线性微分方程用分离变量法可得（2）的通解：y=cep(x)dx（其中c是任意常数） （3）然后从这通解出发，把这通解中的任意常数c编译成的未知函数c(x)，得到y=c(x)e于是：y=c(x)ep(x)dxp(x)dx（4）-c(x)p(x)ep(x)dx（5） p(x)dx将（4）和（5）代入方程（1），得：c(x)ep(x)dx+c(x)p(x)ep(x)dxp(x)dx=p(x)c(x)e+Q(x)即：c(x)e-p(x)dx=Q(x)，所以，c(x)=eQ(x)-p(s)dx所以：c(x)=eQ(x)dx+c所以，（1）的通解为：y=e
11、例1p(x)dx-p(s)dx(eQ(x)dx+c) （其中c是任意常数）dy=-2xy+4x dxdy-x2+2xy=0 解：首先求线性齐次方程的通解y=ce。 dx再应用常数变易法求线性非齐次微分方程的通解，为此，在上式中把常数c变易成待定函数c(x)，即令：y=c(x)e-x，代入原方程得：2c(x)e-x-2xc(x)e-x=-2xc(x)e-x+4x化简得到：c(x)=4xex，上式两边积分得：c(x)=2ex+c22222于是，原方程的通解为y=ce-x+222、一阶非线性微分方程的常数变易法个别的一阶非线性微分方程，可用常数变易法求解，下面介绍四类一阶非线性微分方程的常数变易法。
12、dyyy=+g()x 2.1 齐次方程 dxx（6）对这种方程的解法，在一般教科书中都是首先把它化为可分离变量方程，然后根据可分离变量方程的解法去解，在这里我们可以直接用常数变易法求解。根据常数变易法，先求出原方程“对应”的齐次方程：dyy=的通解：y=cx dxx再令：y=c(x)x （7） 代入（6），有:c(x)x+c(x)=c(x)+gc(x)即：dc(x)gc(x)dc(x)dx=，即： dxxgc(x)x两边积分就可以求出c(x),然后再代入（7），便得原方程的通解。 例2：求方程xy-y=xtan解：将方程改写为yx的通dyyy=+tan dxxxdyy=的通解为：y=cx dx
13、xdc(x)x=tanc(x)，即dx可以求得，它“对应”的齐次线性方程再令：y=c(x)x，代入原方程可得：dc(x)dx=tanc(x)x两边积分得sinc(x)=cx （其中c是任意常数） 代回原变量，得原方程的通解为siny=cx （其中c是任意常数） xdy=p(x)y+Q(x)yn2.2 贝努力方程：dx （8）dy=p(x)y+Q(x)yn形如dx的方程称为伯努利方程,其中p(x),Q(x)为x的连续函数,(n0,1),对于贝努力方程，在一般的教科书上都是先把它化为线性方程，然后根据线性方程的求解方法去解，在这里我们直接用常数变易法去求解。根据常数变易法，先求它“对应”的齐次线性
14、方程令：y=c(x)ep(x)dxp(x)dxdy=p(x)y的通解：y=ce dx代入（8）得，p(x)dxc(x)ep(x)dx+c(x)p(x)en=c(x)p(x)e p(x)dxnp(x)dx+Q(x)cn(x)e即：c(x)=Q(x)c(x)e(n-1)p(x)dx c-n(x)dc(x)=Q(x)e(n-1)p(x)dxdx1n-1(n-1)p(x)dxdx+c解得：c(x)=(1-n)Q(x)e 1n-1所以，（8）的通解为y=ep(x)dx(1-n)Q(x)e(n-1)p(x)dxdx+c 利用此公式可求出任一伯努利方程的通解。 例3、求方程dyy=6-xy2的通解。 dxx
15、6，Q(x)=-x，n=2 x解：可以判断，此方程为贝努力方程，这里p(x)=原方程“对应”的齐次方程为令dyy=6，其通解为：y=cx6， dxxy=c(x)x6，代入原方程化简得：c(x)x6=-xc2(x)x12dc(x)dc(x)=-c2(x)x7，即：-2=x7dx dxc(x)即：1x8=+c c(x)8x6x2所以原方程的通解为：-=c （其中c为任意常数）y8dy=p(x)y2+Q(x)y+R(x)2.3 黎卡提方程：dx （9）一般来说，这一类方程一般来说没有初等解法，不过，若知道其一特解y1，经变换y=z+y1后，方程就变为贝努力方程，因而可解。这里直接用常数变易法求一类特
16、殊的黎卡提方程的解：且a0）-p(x)dxdy，2-2p(x)dx=p(x)y+Q(x)aye+bye+c（a、b、c是实常数，dx根据常数变易法，先求它“对应”的齐次线性方程再令y=c(x)ep(x)dxp(x)dxdy=p(x)y的解y=ce， dx（10）代入原方程，有：dc(x)p(x)dx2e=Q(x)ac(x)+bc(x)+c dx分离变量得到：dc(x)ac(x)+bc(x)+c2-p(x)dx=Q(x)e两边积分，求出c(x)，然后代入（10）可以得原方程的通解。dyy1-x例4、求方程=2+2e(2yex+1)2的通解。dxxx-2dx-p(x)dx1解：在这里由于p(x)=
17、2，e=ex=exx1121ay2e-2p(x)dx+bye-p(x)dx+c=(2yex+1)2 故原方程属于上述黎卡提方程，其中a=4，b=4，c=1。原方程“对应”的齐次线-dyy=2通解为：y=cex 性方程dxx-1x11令y=c(x)e代入原方程有：111211-dydc(x)-x111-=e+c(x)ex(2)=2c(x)ex+2ex(2c(x)exex+1)2 dxdxxxx即：e-1xdc(x)1-x=2e(2c(x)+1)2 dxxdc(x)121-x即： =2edx 22c(x)+1x-1d2c(x)+11x即：=ed- 222c(x)+1x1-1两边积分得：=A-2ex
18、 （其中A是任意常数）2c(x)+11x1x1所以得到：c(x)=e2(Ae-2)所以原方程的通解为：y=1- 212(Ae-2)1x1-ex （其中A为任意常数） 21yy+p(x)e=Q(x)的微分方程 （11） 2.4 形如y先求得（11）“对应”的方程y+p(x)e=0的通解为：y=-lnp(x)dx+c再令：y=-lnp(x)dx+c(x)，代入原方程化简后得： c(x)+Q(x)c(x)=-Q(x)p(x)dx，由此解出c(x)后，便得（11）的通解为-Q(x)dxp(x)dxdx+c y=-lnp(x)dx-eQ(pdx)eyy+p(x)e=Q(x)的通解。 利用此公式可以求得例
19、5、求y+ecosx=y1的通解。 xy-y解：先解方程y+ecosx=0，它的解是e=sinx+c或y=-ln(sinx+c)。可令原方程的解为y=-ln(sinx+c(x)，代入方程得： -c(x)1=sinx+c(x)x即 c(x)+-11c(x)=-sinx xx11xdx1xdxc(x)=e-sinxedx+cx1=-sinxdx+cx1=(cosx+c)xcosxc+) （其中c为任意常数） 所以原方程的通解为y=-ln(sinx+xx 参考文献1、王高雄、周之铭、朱思铭等.常微分方程教程（第三版）M.北京:高等教育出版社.2006.18-642、任永泰、史希福.常微分方程M.沈阳
20、：辽宁人民出版社.1984.4-163、张学元.变系数二阶线性微分方程的一个新的可解类型J.大学数学,2003,19(1):96-98. 4、李鸿祥.两类二阶变系数线性微分方程的求解J.高等数学研究,2002,5(2):10-13. 5、罗亚平,陈仲.微分方程M.南京:南京大学出版社,1987.6、复旦大学数学系.常微分方程M.上海:上海科学技术出版社,1987.7、胡劲松.求一类二阶非齐次欧拉方程的特解J.重庆工学院学报,2003,17(4):49-51. 8、四川大学数学系.高等数学M.北京：人民教育出版社，1979：35-46致 谢在本文的成文过程中，得到了杜老师的悉心指导，在这里我衷心
21、感谢在写论文的这段日子里老师付出的辛苦和所提出的宝贵意见，有了杜老师的不辞辛苦，才有了本文的成形。同时在此，对我所有的专业课任课老师表示崇高的敬意与感谢。另外，由于时间有限，行文较仓促，其中一定还有很多的不足之处，希望各位专家和老师批评指正！谢谢！Ch7-求解常微分方程-1225-232第七章 求解偏微分方程常微分方程的求解是在高等数学里面讲过的，所以大家比较熟悉，对于偏微分方程的求解可能就有点陌生了。事实上偏微分方程在工程上有着更广泛地应用，例如描述液体在多孔介质中的扩散，声学和电磁学中谐波的传播，热在固体中的传导等过程的微分方程都是偏微分方程。所以求解偏微分方程在工程实际上有着非常重要的价
22、值。手工求解常微分方程就已经非常麻烦和复杂了，那么求解偏微分方程就更加的复杂和繁冗了。所幸的是有MATLAB这个计算工具来帮我们解决这一麻烦问题，MATLAB中有一个专门用来求解片微分方程的工具箱PDE工具箱。这个工具箱不但提供有丰富的命令函数，使求解偏微分方程便的简单灵活，而且该提供了一个求解偏微分方程的图形用户界面系统（GUI），使得整个求解过程更加的人性化。特别是对于初学者来说GUIjiang更容易被接受，操作也更方便，所以偏微分方程的求解这一部分我们将主要介绍PDE工具箱的GUI系统。 7-1 偏微分方程的特点对于n阶常微分方程我们知道它的解取决于n个任意常数，例如一阶常微分方程的解可
23、以表示为：u=f(x)+c,c为任意常数。但是对于偏微分方程来说就不是由多少个常数来确2u定的了，例如偏微分方程=f(x,y)的解可表示为： xyu(x,y)=xx0yy0f(x,h)dxdh+w(x)+v(y)其中w(x)和v(y)为两个连续的任意可微函数。我们看到偏微分方程的解可能是非常多的，与常微分方程的解依赖于若干任意常数相比，它的自由度要大的多，对于多维偏微分方程的解更是这样。正是由于这个原因，一般来说偏微分方程的解很难用通式表达出来。事实上我们常用到的是偏微分方程在某种特定条件下的解，这样靠着这些特定条件的约束我们就可以把偏微分方程的解表示出来了。我们把这些帮助确定偏微分方程特解的
24、条件叫做定解条件，由于自变量在多维空间中的变化，其变化的区域非常复杂，所以在区域边界上给出的定解条件也更加形式多样，我们一般称给定在区域边界上的定解条件为边界条件。我们看到边界条件岁于求解偏微由于偏微分方程的求解一般来说太过复杂，所以现在没有一个对所有的偏微分进行求解的理论，所谓的求解偏微分方程也只是对于某些人们比较熟悉的类型的偏微分方程进行求2u2u2u2u解。不同类型的偏微分方程代表不同的物理现象，例如方程=a(2+2+2)表txyz示内部没有热源的情况下物体的温度分布应满足的条件，t表示时间变量，x、y、z分别表示物体的三维空间坐标，这里a=2K，K为物体的热传导系数，Q为该物体的热的容
25、量。Q再比如方程22u2u2u2u=a(2+2+2)表示声波或者振动在介质中的传播，同样t表2txyz示振动传播的时间变量，x、y、z表示介质的三维空间坐标，而a则表示振动在该介质中传播的速度。人们通过对偏微分方程的研究，把常见的偏微分方程归结为四种类型，即：椭圆型方程、抛物型方程和混合型方程。前面讲过的热传导方程是一个抛物型方程，波动方程则是一个双曲型方程。如果热传导方程不考虑时间因素，则方程变为2u2u2ua(2+2+2)=0，这个方程就是一个椭圆方程了，混合型方程顾名思义就是一个xyz2方程具有上面三种类型方程中至少两种的特点的方程了。上面说过，边界条件对于求解偏微分方程很关键，在求解过
26、程中最常用的地定解条件有两种，即：狄利克雷条件（又叫第一类边界条件）、诺曼条件（又叫第二类边界条件）。狄利克雷条件常写为u(t,xi)xW=g(t,xi)，即满足给定区域边界的待求函数u为已知函数g；i对于热传导方程狄利克雷条件的具体意义就是，已知物体的表面温度分布函数为g ;对于波动方程则狄利克雷条件的具体意义就是，波在介质的表面（或介质的给定区域的表面）的振动函数u为已知函数g。诺曼条件常写为un=j(t,xi)，j(t,xi)为已知函数，xiW是n关于的外法线导数。对于热传导方程诺曼条件的具体意义为，沿给定区域边界的外法向的温度变化率为已知函数；对于波动方程则诺曼条件的具体意义为，沿给定
27、区域边界的外法向的振动速率为已知函数。7-2 PDE toolbox求1求解方程的类型使用PDE工具箱求解偏微分方程的基本类型有椭圆型方程、双曲型方程、抛物型方程、特征值方程、椭圆型方程组和非线形椭圆型方程组。PDE工具箱对这些方程的公式表述为：椭圆型方程：-(cu)+au=f,inW; 2u双曲型方程：d2-(cu)+au=f,inW t抛物型方程：du-(cu)+au=f,inW t特征值方程：-(cu)+au=ldu,inW这里是微分算子，表示对多元函数求全微分运算。是给定的平面有界区域，参数c、a、f和待求未知数u都是定义在上的函数。参数d是定义在上的复函数，l是待求的特征值。在抛物型
28、方程和双曲型方程中，参数函数c、a、f、d可以是时间t的函数。可以使用非线形解题器求解非线形椭圆型方程：-(cu)+a(u)u=f(u),这里c、a和f都是未知的待求函数u的函数。另外，PDE工具箱提供的解题器都可以处理下面的方程组：-(c11u1)-(c12u2)+a11u1+a12u2=f1 -(c21u1)-(c22u2)+a21u1+a22u2=f2利用命令还可以求解高阶方程组。2边界条件我们知道常用的两个边界条件狄利克雷条件和诺曼条件，PDE工具箱对它们的公式表示为：Dirichlet条件：hu=r;Neumann条件为：n(cu)+qu=g;这里n是W上的单位外法向矢量，h、r、q
29、和g定义在W上的复值函数（对于特征值问题g=0,r=0）。对于非线性问题，系数g、q、h和 r可以依赖于u；对于抛物型方程和双曲型方程系数还可以依赖于时间t。对于方程组情形，狄利克雷边界条件表示为： rrh11u1+h12u2=r1,h21u1+h22u2=r2一般诺曼边界条件表示为：rrn(c11u1)+n(c12u2)+q11u1+q12u2=g1rrn(c21u1)+n(c22u2)+q21u1+q22u2=g2混合型边界条件表示为：h11u1+h12u2=r1rrn(c11u1)+n(c12u2)+q11u1+q12u2=g1+h11mrrn(c21u1)+n(c22u2)+q21u1
30、+q22u2=g2+h12m这里计算m时要注意满足狄利克雷条件。7-3 GUI求解偏微分方程的过程7-3-1 GUI求解问题的一般程序在命令窗口里输入pdetool，回车，PDE工具箱的图形用户界面（GUI）系统就启动了。 从定义一个偏微分方程问题到完成解偏微分方程的定解，整个过程大致可以分为六个阶段，它们是：（1） 绘制平面有界区域，被MATLAB称之为“Draw模式”。系统提供了一系列的实体模型，包括矩形、圆、椭圆和多边形，我们可以通过公式把这些实体模型组合起来，生成我们需要的平面区域。（2） 定义边界，被MATLAB称之为“Boundary模式”。声明不同边界段的边界条件。（3） 定义偏
31、微分方程，被MATLAB称之为“PDE模式”。确定方程的类型和方程的系数c、a、f和d。根据具体情况（例如材料的性质）,我们还可以在不同的子区域上独立的声明不同的系数。（4） 网格化区域，被MATLAB称之为“Mesh模式”。我们可以控制自动生成网格的参数，还可以对生成的网格进行对次的细化，使网格分割的更细更合理。（5） 解偏微分方程，被MATLAB称之为“Solve模式”。对于椭圆型方程，我们可以激活并控制非线性自适应解题器来处理非线性方程；对于抛物线方程和双曲线方程，设置初始边界条件后可以求出给定时刻t的解；对于特征值问题，可以求出给定区间上的特征值。求解完成后我们还可以返回到第四步，对网
32、络进一步细化，进行再次求解。（6）计算结果的可视化。我我们可以通过设置系统提供的对话框，显示所求的解的表面图、网络图、等高线图和箭头梯度图（quiver）。对于抛物型和双曲型问题还可以进行动画演示。如何通过上面的六个阶段来具体的求解一个偏微分方程，我们将在下一节中通过一个示例给出具体的操作步骤。 7-3-2 求解问题的示例让我们从求解最简单的偏微分方程泊松方程开始，泊松方程的表达式为：2u2u2u-Du=-(2+2+2)=f xyz这里D是拉普拉斯算子，注意和微分算子区分开。定解条件我们采用狄利克雷边界条件。我们看到泊松方程属于椭圆型偏微分方程，它可以表示内部没有热源的物体的热的分布，还可以表
33、示静止的电磁场的分布。首先在MATLAB的命令窗口输入pdetool，回车确认，这样图形用户界面系统（GUI）就启动了。接下来我们将按照上一节所说的六个步骤开始操作了，启动后的GUI窗口如图7-1所示。图7-1 PDE工具箱的GUI 窗口 在打开的窗口中，点击菜单Options，选中Grid选项，于是窗口中便出现了网格，这样方便我们在绘图时确定几何图形的位置和大小。用鼠标左键单击工具栏里的画圆按钮（左边第四个），使其处于选中状态（表现为下沉），然后在窗口的绘图区域按下左键拖动鼠标，于是一个椭圆就出现了，系统给它命名为E1。再选中工具栏里的画矩形的按钮（左边第二个），在绘图区域的中心位置拖动鼠标
34、，于是一个矩形就出现了，系统给它命名为R1。把光标移到E1上面，双击鼠标左键，这时会弹出一个对话框，用来定义E1的位置和大小，设置如图7-2所示。 图7-2 几何实体E1属性对话框 然后选择OK按钮确认。用同样的办法调出R1的属性对话框，设置如下，然后确认，如图7-3所示。 图7-3 R1属性对话框 接下来把窗口中Set formula栏中的公式修改为E1-R1，这样几何实体的绘制就完成了，画完以后窗口如图7-4所示。我们用鼠标点住画好的几何实体，然后拖动鼠标就可以移动实体的位置，另外还可以把画好的几何实体用m文件的形式保存起来。图7-4 绘好的几何实体接下来设置边界。用鼠标左键单击工具栏里的
35、按钮W，这样边界就生成了，如图7-5所示。 图7-5 生成边界条件我们看到边界是分段的并且有方向（带有箭头）。把光标移到一段边界上，然后双击鼠标左键，就会出现一个对话框，让我们设置该段的边界条件，如图7-6所示。本例中所有的边界都采用Dirichlet条件的默认值。 图7-6 边界条件对话框我们既可以采用上面的办法一段一段的设置边界条件，也可以一次设置好边界条件相同的边界段。先按下shift键，然后用鼠标左键分别单击条件性同的边段，这样这些边界段都被选中了，再使用上面的办法弹出对话框设置参数，于是被选中的段都被设置了。选择方程的类型。单击工具栏里的PDE 按钮，打开方程 类型对话框，由于 泊松
36、方程适宜中椭圆型方程，所以在对话框中选择Elliptic选项。方程的参数设置如图7-7所示，单击OK按钮确认。 图7-7 方程类型对话框对设定的区域进行网格剖分。单击工具栏里的三角号按钮，就可以进行网格剖分了。我们还可以单击工具栏里网格剖分按钮右边的按钮（一个三角号里面还有一个小三角的按钮）对网格剖分进一步细化，网格细化后的GUI窗口如图7-8所示。图7-8 网格细化后的图像求解方程。单击工具栏里的等号按钮，开始求解。这时窗口里就显示出了泊松方程在给定边界条件下的图形解，如图7-9所示。图7-9 泊松方程的图形解该图中的图像是通过颜色的深浅来表示z坐标的大小的。 我们不但能用平面图来表示方程的
37、解，还可以使用多种方式可视化方程的解。例如在图7-9中，用鼠标单击等号按钮右边的那个图标，就会弹出一个对话框，在这个对话框里可以选择我们想要的解的图像表现形式。设置好该对话框，然后单击plot按钮，我们想要的结果就出现了。本例的对话框设置和给出的图像如图7-10、7-11所示。图7-10 绘图对话框在图7-11中，我们可以用鼠标点住图像然后拖动它，从不同的角度观察解的情况。 6常微分方程的求解6常微分方程的求解一、知识背景常微分方程在物理学，工程技术中运用非常广泛，相当重要。许多物质运动的过程用常微分方程来描述，如：质点的加速运动、谐振动、电容充放电过程及电感通电断电等过程，因此，求解常微分方
38、程成为很多物理问题求解的一种常用方法。而有时很难求出解析解，但可求出常微分方程的数值解以逼近解析解，以完成对摸型的研究。求解常微分方程数值解通常有欧拉法和龙格-贝塔法。欧拉法解法的基本思想是在小区间上用差商代替导数，并通过把小区间不断地划分求极限，从而最终得到数值解。龙格-贝塔方法基本思想也是用差商代替导数，但是该方法是在小区间内运用了微分中值定理，在小区间内多取点，再取加权平均值来构造精度更高的计算式。在MATLAB系统中，主要采用龙格-贝塔法来计算常微分方程的数值解。 举例：如图6-1中所示电路，先将开关k掷向“1”端，待电容器c充完电后，将开关k掷向“2”端，电容器开始放电。放电过程满足
39、下面方程：Rdqdt+qc=0 图 6-1dqq=- dtcdqdt=f(q,t)或 R我们也可写成右面方程组形式：，q(t0)=q0这即是一阶常微分方程初值问题的一般形式。二、 计算指令 ode23，ode45语句格式（以ode23为例）： t, y = ode23 (f,tspan,y0,tol )语句中各符号意义如下f：求解的常微分方程的文件名，把方程写成函数形式并存储于m文件中。方程形式为y=f(t,y)。 115tspan：输入t0,tf,分别为自变量的初始值和最终值，为单调递增（减）的积分区间。y0：函数的初始值矢量。Tol：误差范围，(缺省值为0.000001)t,f：t是输出的
40、时间列矢量，矩阵y的每个列矢量是解的一个分量。 例1：用求数值解方法,求解常微分方程：x=-x3，初始值x(0)=1。 解：定义并储存函数文件: yfun.m function exer=yfun(t,x);exer=-x.3; 程序文件：t,x=ode23(yfun,0,1,1);plot(t,x,k+,t,x,-r)例 1 图* 请改指令ode45,看图形有什么区别。Ode23，低精度，使用Runge-Kutta法的二，三阶算法。 Ode45，中等精度，使用Runge-Kutta法的四，五阶算法。三、求解二阶以上的常微分方程的数值解将方程在函数中迭代表达，化成一阶微分方程组，从而在命令行中
41、编程实现。下面以单摆模型为例，说明二阶以上常微分方程的一般解法。 例2：如下图所示，在单摆模型中，以铅直方向为平衡位置，以右边为正方向建立摆角x的坐标系。 .x0=0.1745(rad), w=x=0(rad/s)。116解：依据转动定律vvvvM=rF=Jbmld2dxdt22=-lmg(sinx)glsinx 2x2dt=- MATLAB令：y (1) = x，y (2)=dx/dt，则：dy(1)dtdy(2)dt=y(2)g=sinxl 图 6-2函数文件：（保存yfun2）function exer2=yfun2(t,y);g=9.8;length=25;exer2=y(2);-g/
42、length*sin(y(1); 程序文件：ts=0,10;y0=0.1745,0;t,y=ode45(yfun2,ts,y0); subplot(1,2,1), plot(t,y),grid, subplot(2,2,2),plot(y(:,1),y(:,2),grid例 2 图 * 关于函数文件的讨论： 117必须严格按照格式编写，不得随意变动！在函数文件中，要求两个等号左边的变量名相同，在此用的是exer2。对比原方程知：exer2(1)=dy(1)/dt，exer2(2)=dy(2)/dt。等号右边yfun2为存盘文件名，在程序文件中，必须以这名称提取。在括号中，t为时间变量，y为在一
43、阶常微分方程组定义的量。 * 关于程序文件的讨论：ode45为解常微分方程组的指令。 调用文件，必须用函数名yfun2调用。 时间范围从0刭10，步长为默认值。 初始角位移为0.1745，初始角速度为0。输出结果为时间变量的值t、角位移的值y(:,1),角速度的值y(:,2)（需把语句中；号改成，号才可显示）。最后用plot指令绘制函数曲线y(1)t（即xt），相图y(2)y(1)图线。 可增添语句，了解更多的东西，（请作为练习）如，键入size(t),size(y)-可知t,y的大小，各矢量的长度，即元素个数。 键入t,y，可得t,y的数值解。 * 关于物理问题的一些讨论（1）一条曲线是xt
44、，另一条是xt，后者是前者的时间变化率。 （2）可以在同一图中作一条余弦曲线，（请试作一下）再和xt比较，改变角位移初始值，从不同曲线了解x0的影响。（3）关于相图，角位移 和角速度 是两个动力学变量，它们构成一个相平面，相平面中的每一点代表系统的一种可能的运动状态。封闭轨道相当于周期振动。（4）再看一下阻尼摆：djdt22+kdj+wdt2sin(j)=0例3：同样，令：y（1）= x，y（2）= dx / dtdy(1)=y(2)dtdy(2)=-gsin(y(1)-0.1*y(2)ldt（设k=0.1）function xexer5=yfun5(t,y);g=9.8;l=25;xexer
45、5=y(2);-g/l*sin(y(1)-.1.*y(2); 118ts=0,10;y0=0.1745,0;t,y=ode23(yfun5,ts,y0);subplot(1,2,1), plot(t,y),grid,subplot(1,2,2),plot(y(:,1),y(:,2),grid 从上图可见阻尼摆在相平面中的运动轨道不再是周期的，而是盘旋缩小，最终静止在原点的一条螺旋线。例4：设t = 0100s，再绘制阻尼摆的图线。(下图，ts=0,100)例 3 图例 4 图 例5：为了看到混沌运动，可以再对摆施加周期性外力，使它成为强迫摆。 119djdt22+kdjdt+w2(1+Asin
46、(Wt)sinj=0，设定：k=0.1，A=1，W=1。function：xexer5=yfun5(t,y);g=9.8;l=25;xexer5=y(2);-g/l*sin(y(1).*(1+sin(t)-.1.*y(2); ts=0,95;y0=0.1745,0;t,y=ode23(yfun6,ts,y0);subplot(1,2,1),plot(t,y),grid,subplot(1,2,2),plot(y(:,1),y(:,2),grid例 5 图 四、小结求解常微分方程主要分三步：编写函数文件（odefile）并储存；选择恰当的解方程的指令，指定解方程的耍求；编写程序按规定的格式调用解方程组。在MATLAB6.0中有关编写odefile模板，可以查看。 more on, type odefile, more off求解方程组的指令都有默认的求解的条件和要求，可以用指令odeset修改或重新确定求解的条件和要求。可以用指令odegest获取优化选项中参数取值信息。 120解方程组指令更普遍语句格式： t, f = ode23 (f, tspan