mathematica.doc

第二讲 微分方程实验一、实验目的现实世界的运动规律，很多都是用微分方程模型来描述，大到天体运动，小到基本粒子的运动，无一不是如此.因此、微分方程在力学、物理学、各种工程技术科学、生物生态学以及经济学等领域都有非常重要的应用.这里涉及两类问题，一类是如何将实际问题转化为微分方程模型，另一类是如何求出方程的解.本实验将借助Mathematica研究与一阶常微分方程相关的几个问题，让我们了解几类微分方程模型、利用Mathematica求解微分方程的方法及如何直观演示方程解的几何形态.二、基本理论让我们首先复习一阶微分方程（组）的几个概念以及利用Mathematica求微分方程解的方法.对于一阶微分方程， （1）若存在定义在某区间 内的可微函数 满足，则称 为方程（1）的解， 所对应的曲线称为方程的解曲线（或积分曲线），它直观反映了方程解的变化.若解 满足条件（2）则称 为满足初始条件（2）的特解，此时解曲线经过点 .同样，对于一阶微分方程组， （3）若存在定义在某区间 内的可微函数 满足则称 为方程组（2）的解，若还满足， （4）则称 为满足初始条件（4）的特解.曲线（5）称为方程组（3）解的相轨线，而解曲线是三维空间 中的曲线（6）若方程（3）在无穷区间 内有解，且满足解的存在唯一性条件，则（3）的解对应一个动力系统.研究动力系统在无穷远处的变化性态是非常有意义的事情.若方程（1）的解 或方程组（3）的解 能用初等函数表示，则称其为各自方程的解析解.我们在微分方程的求解理论中了解了求方程解的一些方法，利用Mathematica 也可以完成求解过程.首先让我们回顾一下Mathematica中用于求微分方程通解和特解的函数及其用法.DSolveyxfx,yx,yx,x-求方程（1）的通解DSolveyxfx,yx,yx0y0,yx,x-求方程（1）在初始条件（3）下的特解DSolvextft,xt,yt,ytgt,xt,yt,xt,yt,t-求方程组（4）的通解DSolvextft,xt,yt,ytgt,xt,yt,t0x0,yt0y0,xt,yt,t-求方程组（4）在初始条件（6）下的特解NDSolveyx=fx,yx,yx0y0,yx,x,x0,x1-求方程（1）在初始条件（3）下数值解 NDSolvextft,xt,yt,ytgt,xt,yt,xt0x0,yt0y0,xt,yt,t,t0,t1-求方程组（4）在初始条件（6）下的数值解特别提示 在微分方程求解的操作中，一定要将未知函数y写成yx或yt的形式！另外，用NDSolve求方程的数值解时，其输出形式是插值函数InterpolatingFunction.例1 利用Mathematica求微分方程 的通解及在初始条件 下的特解.MathematicaDSolveyx=(1+x2) yx,yx,xDSolveyx=(1+x2) yx,y0=1,yx,x运算结果分别为：在实际问题中，我们往往要利用方程的解，如画出解曲线的图形，那么如何将上述结果变为可用的形式？我们可以利用替换的方式将方程解的表达式转化为表的形式：yx/.DSolveyx=(1+x2) yx,yx,xyx/.DSolveyx=(1+x2) yx,y0=1,yx,x运算结果为：此时可以通过对表的操作引用方程的解.例2 利用Mathematica求微分方程组 通解和在初始条件 下的特解.DSolvext3xt-2yt,yt2xt-yt,xt,yt,tDSolvext3xt-2yt,yt2xt-yt, x01,y00,xt,yt,t运算结果为：但多数情况下没有这么幸运，经常碰到的情形是方程不存在解析解，或者我们能人工求解Mathematica却无能为力.读者不妨尝试一下用上面的方法求黎卡提（Riccati）方程和方程组的解，观察会发生怎样的情形.你会看到第一个方程的一个形式非常复杂的解，实际上，它是用特殊函数来表示的，而非通常的初等函数，因此已不是方程的解析解.而对于第二个方程组，正确地输入以后，运行所得到的结果与输入同样的形式，什么也没有做！也许你对Mahematica的能力产生了怀疑，实际上我们即使不能求得解析解，借助它来研究方程解的性质.在这一章中，我们将介绍如何描绘微分方程所确定的线素场，如何得到方程的几何和数值近似解，同时给出微分方程在实际问题中的几个有趣的应用.三、实验内容问题一：磁场的分布-微分方程的线素场如果你手中有一个条形磁铁和少许铁质粉末，你便可以通过实验了解在磁铁周围的平面内的磁场分布情况，下面我们讨论与此相关的问题.微分方程 告诉了我们这样的信息，在平面上某区域 内一点 处，解曲线 的斜率为 .我们在区域 内的有限个点处作从该点出发的短线段，使其斜率为函数 在该点的值，这些线段构成的集合称为微分方程的线素场（line element field）或方向场（directional field），由于每一线段均与方程的解曲线相切，因此当线段长度很短时，它们是解曲线的局部近似，所以线素场可以反映解的变化趋势，当所选取的点俞密，这种变化趋势俞明显.那么，如何用Mathematica来绘制方程的线素场呢？让我们来寻找解决问题的方法.首先我们作出从点 出发，斜率为 且长度为 的线段.我们只需在过 且斜率为 的直线上找到一点 ,使 与 的距离为 ,利用 Mathematica 解方程且取其中一组解即可.Mathematica程序（ch2-ex1.nb）运算结果为：因此利用基本图形Line便可做出所需线段，例如给定 ，即a=1;b=2;h=1;k=1/2;GraphicsLinea,b,xa,b,h,k,ya,b,h,k/Show输出结果为：下面我们生成微分方程的线素场.首先取函数 定义域中的闭矩形区域，将区间 和 分别 等分和 等分，得区域 的一个划分，即区域 的网格图,将网格的每个交点作为小线段的出发点，这样便得到方程分布在区域 中的线素场.我们结合前面画线段的方法，可以编写如下Mathematica程序模块.Mathematica程序（ch2-ex2.nb）ClearAllf,hdirectionalFieldf_,x_,xmin_,xmax_,y_,ymin_,ymax_,m_,n_,h_,opts_:=Modules,s1,x,y,a,b,k,i,j,lines,s=Solvey-b=k (x-a),(x-a)2+(y-b)2=h2,x,y/Simplify;s1=Flattens;xa_,b_,k_=x/.s14;ya_,b_,k_=y/.s13;h1=(xmax-xmin)/m;h2=(ymax-ymin)/n;a=xmin+i h1;b=ymin+j h2;k=f/.x-a,y-b;lines=TableLinea,b,xa,b,k,ya,b,k,i,0,m,j,0,n;ShowGraphicslines,opts,Axes-Automatic,AspectRatio-Automatic通过上述程序模块，我们可以对具体的方程描绘它的线素场，让我们看两个例子.fx_,y_:=x2+y2-1directionalFieldfx,y,x,-2,2,y,-2,2,20,20,0.1运行结果：fx_,y_:=3/2Signy Absy(1/3)directionalFieldfx,y,x,-2,2,y,-2,2,20,20,0.1运行结果：最后，我们观察一个实际问题的线素场，可以看到其几何描述的正好与实际情形相吻合. 在平面 和 处分别放置两个正、负单位电荷，则它们在平面上产生一磁场. 从微分方程教材中，获知描述磁场强度的微分方程为其中对具体的 ，我们可以作出磁场的分布图.Mathematica程序（ch2-ex3.nb）a=2;Px_,y_:=(x+a)/(x+a)2+y2)(3/2)-(x-a)/(x-a)2+y2)(3/2)Qx_,y_:=y/(x+a)2+y2)(3/2)-y/(x-a)2+y2)(3/2)fx_,y_:=Qx,y/Px,ydirectionalFieldfx,y,x,-3,3,y,-3,3,20,20,0.15,GraphicsRGBColor1,0,0,PointSize0.03,Point-2,0,Point2,0运行结果：这里我们看到在程序模块directionalField中定义变量opt给我们带来了方便，通过它我们可以丰富线素场图，如上面画出两个点电荷的位置，使得到的图形更形象，读者可以在有关的教材中寻找更多的例子.问题二：炮弹飞行的轨迹-欧拉方法在炮弹发射时，我们如何知道炮弹在空中飞行中的轨迹? 假设你拥有一架高速照相机，你可以将炮弹飞行的过程记录下来.但冲洗出来的照片并不是一条连续的轨线，而是炮弹在不同位置的离散图象，不过这已足够让我们了解炮弹的飞行轨迹了，只要照相机的速度有足够的快.前面我们已经知道很多情况下我们不能求出微分方程的解析解，但如照相一样，我们只要在解曲线上找到一连串的点，也可以获知解曲线的形状. 在前面的实验中我们了解了描绘一阶常微分方程线素场的方法，由线素场我们可大概了解解曲线的走势，但它并不能较准确地反映每一条解曲线的几何形状,如何在不能求出方程的解析解的情况下，能将解曲线比较精确地描绘出来？欧拉为我们提供了一个方法，其原理是先通过差分将一阶常微分方程化为差分方程（即代数方程），然后考虑代数方程的解（点列），将这些点顺次连接起来，变得到解曲线的近似图形.考虑初值问题由导数的定义知，当 时， ，因此上述方程可改写成，或 取一列点 ，记 ，则有，我们称该方程为差分方程， 称为差分步长，一般来说， 越小，所求得的点列越反映微分方程解的实际情况，这种处理微分方程的方法称为欧拉（Euler）方法.首先，我们给出用欧拉方法求微分方程近似解的Mathematica程序.Mathematica程序（ch2-ex4.nb）ClearAllfeulerf_,x_,x0_,x1_,h_,y_,y0_,opt_:=Modulepoints,y0=y0;xn_:=x0+n h;yn_:=yn-1+ h f/.x-xn-1,y-yn-1/N;m=Floor(x1-x0)/h;points=Tablexn,yn,n,0,m;ListPlotpoints,opt我们取 ，运行下面的程序fx_,y_:=3/2 Signy Absy(1/3)eulerfx,y,x,0,1,0.1,y,1/2,AxesOrigin-0,0,PlotRange-2.2,2.2,-2.2,2.2,PlotStyle-RGBColor1,0,0,PointSize0.02,AspectRatio-Automatic运行结果：结合实验一中的线素场，我们可以看到欧拉方法所得解与线素场的一致性.Mathematica程序（ch2-ex5.nb）f1=eulerfx,y,x,0,1,0.1,y,1/2,AxesOrigin-0,0,PlotRange-2.2,2.2,-2.2,2.2,PlotStyle-RGBColor1,0,0,PointSize0.02,AspectRatio-Automatic,DisplayFunction-Identity;f2=directionalFieldfx,y,x,-2,2,y,-2,2,20,20,0.1,DisplayFunction-IdentityShowf1,f2,DisplayFunction-$DisplayFunction运行结果：另外，我们可以将方程的解析解与欧拉方法所得到的近似解进行比较，让我们以初值问题为例作出解析解和近似解的图形.Mathematica程序（ch2-ex6.nb）Clearyfx_,y_=Sinx+y;s=DSolveyx=fx,yx,y0=1,yx,x;p1=Plotyx/.s1,x,0,2,DisplayFunction-Identity;p2=eulerfx,y,x,0,2,0.2,y,1,PlotStyle-RGBColor1,0,0,PointSize0.02,DisplayFunction-Identity;Showp1,p2,DisplayFunction-$DisplayFunction;运行结果：上面介绍的欧拉方法可以改进，使得到的结果更精确，基本想法是取两个斜率的平均值，即我们通过具体例子比较改进前后的欧拉方法的计算结果和精确解的结果，请读者考虑初值问题编写出获得上述结果的Mathematica程序.问题三：黄土高原上的沟壑-最陡下降法也许你去过西北的黄土高原，也许从各种媒体上了解那里的景象，到处是沟壑交错，这主要是由于长年雨水的冲刷而形成的. 你也许会想到，这些沟壑的分布走向与地形有关，如果你手中有一张当地的地形图，你会发现这些沟壑往往是垂直穿过所倚山峦的等值线，下面我们将用数学的理论给出解释.假设山坡可以用函数 来表示，且坡面是光滑的，若落在山坡上的雨水没有渗透，则它应往最陡的方向往下流动，现在我们确定流动的路线.我们考虑路线在山脚所在平面（即 平面）投影，设投影曲线方程为 ,则由方向导数和梯度的意义，可建立方程（7）以该投影曲线为准线作母线垂直于z轴的柱面，则该柱面与坡面的交线便是雨水在坡面流下的路线.下面我们给出计算实例.设山坡是上半椭球面,雨水落下的位置为 ,我们利用Mathematica完成上面的任务,分下面几步进行：（1） 求投影曲线方程.设 ， 为为其在平面的投影曲线方程，解初值问题（7）便求得投影曲线方程.不妨设 在第一卦限.（2） 求出投影曲线的自变量的变化范围.解方程组，解得正解 ，则 得变化范围为 .（3） 描绘雨水下流路线及投影曲线，并从不同的角度观察图形.Mathematica程序（ch2-ex7.nb）quickestdescenta_, b_, c_, x0_, y0_, h1_, h2_, opt_:=Modulef, fx, fy, g, ellipsoid, s, t, r, l, p, d, descentroad, projection, generator, y,fx_, y_:=c Sqrt1-x2/a2-y2/b2;fxx_, y_:=Dfx, y, x;fyx_, y_:=Dfx, y, y;gx_, y_=fyx, y/fxx, y;ellipsoid=ParametricPlot3Da Sinu Cosv, b Sinu Sinv, c Cosu,u, 0, Pi/2, v, 0, 2Pi, DisplayFunction-Identity;s=DSolveyx=gx,yx, yx0=y0, yx, x;yx=yx/. s1,1;Printy(x)= , yx;t=Solvefx,y=0, y=yx/. s1,1, x, y;Fori=1, i0, l=r;p=Tablex, yx/. s1,1, 0, x, x0, l, h1;projection=Graphics3DRGBColor1,0,0, Thickness0.015, Linep;d=Tablex, yx/. s1,1, fx, yx/. s1,1, x, x0, l, h1;descentroad=Graphics3DRGBColor1,0,0,Thickness0.015, Lined;generator=TableGraphics3DRGBColor0, 1, 0,Linex, yx/. s1,1, 0, x, yx/.s1,1, fx, yx/.s1,1, x, x0, l, h2;startpoint=Graphics3DPointSize0.03, RGBColor1,0,0,Pointx0, y0, fx0, y0;Showprojection, ellipsoid, startpoint, descentroad, generator,PlotRange-a,a,0,b,0,c,opt给出具体参数，便可得到非常直观得结果.例如x0=0.1;y0=1.5;a=3;b=6;c=3;f1=quickestdescenta,b,c,x0,y0,0.001,0.2,ViewPoint-1,-2,1,DisplayFunction-Identity,Boxed-Falsef2=quickestdescenta,b,c,x0,y0,0.001,0.2,ViewPoint-1,2,1,DisplayFunction-Identity,Boxed-Falsef3=quickestdescenta,b,c,x0,y0,0.001,0.2,ViewPoint-1,-2,0.7,DisplayFunction-Identity,Boxed-FalseShowGraphicsArrayf1,f2,f3运行结果：问题四：战争能爆发吗？-微分方程的稳定性两国势均力敌的军事力量互相制约，能保证相互之间的和平共处，一旦某一国的军事力量无限扩张，便会对另一国构成威胁，爆发战争的机会大大增加.在建立两国军备竞争的数学模型之前，让我们观察几类微分方程的渐近性态，掌握了解方程变化趋势的几何方法.例1 考虑方程 的解在平衡点 附近的性态.Mathematica程序（ch2-ex8.nb）Clearx,y,s,datas=TableNDSolvext=-xt-yt,yt=xt-3yt,x0=x0,y0=0,x,y,t,0,20,x0,-1,1,0.2;a=x,y/.s;xi_:=ai,1,1yi_:=ai,1,2orbit=TableParametricPlotxit,yit,t,0,20,PlotRange-0.8,0.8,-0.2,0.2,PlotStyle-RGBColor1,0,0,Thickness0.01,DisplayFunction-Identity,i,1,10;Showorbit,DisplayFunction-$DisplayFunction;运行结果：此时，在不同初值下的轨线趋于平衡点，称平衡点是稳定的.我们还可以通过下面的数值得到相应的结果.Mathematica程序（ch2-ex9.nb）Clearx,y,s,datas=TableNDSolvext=-xt-yt,yt=xt-3yt,x0=x0,y0=0,x,y,t,0,20,x0,-1,1,0.2;a=x,y/.s;xi_:=ai,1,1yi_:=ai,1,2data=Tablet,x1t,y1t,t,0,20,1;GridBoxJoint,x(t),x(0)=-1,y(t),y(0)=0,data, ColumnLines-True,RowLines-True,ColumnSpacings-3,5/DisplayForm运行结果：从上述数据表格中我们发现，随着 的增加， 和 均趋于零,这与前面的图形正好吻合.是否这一类方程的解都具有如此性质，我们不妨再看一个例子.例2考虑方程 的解在平衡点 附近的性态.Mathematica程序（ch2-ex10.nb）Clearx,y,s,datas=TableNDSolvext=xt-yt,yt=2xt+yt,x0=a,y0=0,x,y,t,0,20,a,-1,1,0.2;a=x,y/.s;xi_:=ai,1,1yi_:=ai,1,2data=TableParametricPlotxit,yit,t,0,20,PlotRange-20,20,-20,20,PlotStyle-RGBColor1,0,0,Thickness0.01,DisplayFunction-Identity,i,1,10;Showdata,DisplayFunction-$DisplayFunction;运行结果：此时当 增大时， 将无限远离平衡点，这时我们称原点是非稳定的平衡点，同样也可以用数据加以说明.Mathematica程序（ch2-ex11.nb）Clearx,y,s,datas=TableNDSolvext=xt-yt,yt=2xt+yt,x0=a,y0=0,x,y,t,0,20,a,-1,1,0.2;a=x,y/.s;xi_:=ai,1,1yi_:=ai,1,2data=Tablet,x1t,y1t,t,0,20,1;GridBoxJoint,x(t),x(0)=-1,y(t),y(0)=0,data,ColumnLines-True,RowLines-True,ColumnSpacings-3,5/DisplayForm运行结果：现在我们回到开始的问题.一个国家的军事力量与该国的军队数量、武器装备和经济实力等因素有关，我们将甲乙两国的军事力量量化为 和 .一般来说，影响军事力量增加的速度有三个方面：