1.画出下列微分方程的向量场和一些积分曲线.

Mathematica 与常微分方程—方向场和积分曲线摘要：长期以来，从小学到大学十几年，数学一直是我们学习的一门主课，老师所讲的、学生所练、所考的主要是定义叙述、定理证明、公式推算、计算方法、……，数学给我们的印象是，沿定义→公理→定理→推论→证明这么一条演绎道路进行的、一个十分严格的数学推理王国和一个充满美感的抽象世界。

然而，我们却不知道，也许也没有想过，这些如此严密、完整、美妙的结论是怎么来的？数学家是通过什么样的方式发现它们的？我们从这些可爱结论本身看不到数学家发现它们的艰辛，也体会不到数学家在发现它们之后的一种喜悦。

关键词：Mathematica常微分方程方向场积分曲线简介Mathematica 是美国Wolfram Research公司研制的一种数学软件，集文本编辑、符号计算、数值计算、逻辑分析、图形、动画、声音于一体，与Matlab、Maple 一起被称为目前国际上最流行的三大数学软件。

它以符号计算见长，同时具有强大的图形功能和高精度的数值计算功能。

在Mathematica中可以进行各种符号和数值运算，包括微积分、线性代数、概率论和数理统计等数学各个分支中公式的推演、数值求解非线性方程、最优化问题等，可以绘制各种复杂的二维和三维图形，并能产生动画和声音。

Mathematic系统与常见的高级程序设计语言相似，都是通过大量的函数和命令来实现其功能的。

要灵活使用Mathematica，就必须尽可能熟悉各种内部函数（包括内置函数和软件包函数）。

由于篇幅限制，本附录仅分类介绍Mathematica的基本功能，及与微积分有关的函数（命令）的使用，其他功能请读者自行查阅帮助或有关参考文献。

方程对于学过中学数学的人来说是比较熟悉的；在初等数学中就有各种各样的方程，比如线性方程、二次方程、高次方程、指数方程、对数方程、三角方程和方程组等等。

这些方程都是要把研究的问题中的已知数和未知数之间的关系找出来，列出包含一个未知数或几个未知数的一个或者多个方程式，然后取求方程的解。

实验2--微分方程(基础实验)119 项目四 无穷级数与微分方程实验2 微分方程（基础实验）实验目的 理解常微分方程解的概念以及积分曲线和方向场的概念,掌握利用Mathematica 求微分方程及方程组解的常用命令和方法.基本命令1. 求微分方程的解的命令DSolve对于可以用积分方法求解的微分方程和微分方程组,可用Dsolve 命令来求其通解或特解.例如,求方程023=+'+''y y y 的通解, 输入DSolve[y ''[x]+3y '[x]+2y[x]==0,y[x],x]则输出含有两个任意常数C[1]和C[2]的通解:{}{}]2[C e ]1[C e ]x [y x x 2--+→注:在上述命令中,一阶导数符号 ' 是通过键盘上的单引号 ' 输入的,二阶导数符号 '' 要输入两个单引号,而不能输入一个双引号.又如,求解微分方程的初值问题:,10,6,03400='==+'+''==x x y y y y y输入Dsolve[{y''[x]+4 y'[x]+3y[x]==0,y[0]==6, y'[0]==10},y[x],x](*大括号把方程和初始条件放在一起*)则输出{}{}x 2x 3e 148(e ]x [y +-→-2. 求微分方程的数值解的命令NDSolve对于不可以用积分方法求解的微分方程初值问题,可以用NDSolve 命令来求其特解.例如要求方程5.0,032=+='=x y x y y的近似解)5.10(≤≤x , 输入NDSolve[{y'[x]==y[x]^2+x^3,y[0]==0.5},y[x],{x,0,1.5}](*命令中的{x,0,1.5}表示相应的区间*)则输出{{y->InterpolatingFunction[{{0.,1.5}},< >]}}注:因为NDSolve 命令得到的输出是解)(x y y =的近似值. 首先在区间[0,1.5]内插入一系 列点n x x x ,,,21Λ, 计算出在这些点上函数的近似值n y y y ,,,21Λ, 再通过插值方法得到 )(x y y =在区间上的近似解.3. 一阶微分方程的方向场一般地,我们可把一阶微分方程写为),(y x f y ='的形式,其中),(y x f 是已知函数. 上述微分方程表明:未知函数y 在点x 处的斜率等于函数120f 在点),(y x 处的函数值. 因此,可在Oxy 平面上的每一点, 作出过该点的以),(y x f 为斜率 的一条很短的直线(即是未知函数y 的切线). 这样得到的一个图形就是微分方程),(y x f y ='的方向场. 为了便于观察, 实际上只要在Oxy 平面上取适当多的点,作出在这些点的函数的 切线. 顺着斜率的走向画出符合初始条件的解,就可以得到方程),(y x f y ='的近似的积分曲 线.例如, 画出0)0(,12=-=y y dxdy 的方向场. 输入<<Graphics`PlotField`g1=PlotVectorField[{1,1-y^2},{x,-3,3},{y,-2,2}, Frame->True,ScaleFunction->(1&),ScaleFactor->0.16,HeadLength->0.01,PlotPoints->{20,25}];则输出方向场的图形（图2.1）, 从图中可以观察到, 当初始条件为2/10=y 时, 这个微分方程的解介于1-和1之间, 且当x 趋向于-∞或∞时, )(x y 分别趋向于1-与1.-3-2-10123-2-1012 -3-2-10123-2-112下面求解这个微分方程, 并在同一坐标系中画出方程的解与方向场的图解. 输入sol=DSolve[{y'[x]==1-y[x]^2,y[0]==0},y[x],x];g2=Plot[sol[[1,1,2]],{x,-3,3},PlotStyle->{Hue[0.1],Thickness[0.005]}];Show[g2,g1,Axes->None,Frame->True];则输出微分方程的解xxe e x y 2211)(++-=,以及解曲线与方向场的图形（图2.2）. 从图中可以看到, 微分方程的解与方向场的箭头方向相吻合.实验内容用Dsolve 命令求解微分方程例2.1 (教材 例2.1) 求微分方程 22x xe xy y -=+'的通解.输入Clear[x,y];DSolve[y '[x]+2x*y[x]==x*Exp[-x^2],y[x],x]或DSolve[D[y[x],x]+2x*y[x]==x*Exp[-x^2],y[x],x]则输出微分方程的通解:121 ⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧+→--]1[C e x e 21]x [y 22x 2x 其中C[1]是任意常数.例2.2 (教材 例2.2) 求微分方程0=-+'x e y y x 在初始条件e y x 21==下的特解. 输入Clear[x,y];DSolve[{x*y ' [x]+y[x]-Exp[x]==0,y[1]==2 E},y[x],x]则输出所求特解:⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧+→x e e ]x [y x 例2.3 (教材 例2.3) 求微分方程x e y y y x 2cos 52=+'-''的通解.输入DSolve[y ''[x]-2y '[x]+5y[x]==Exp[x]*Cos[2 x],y[x],x]//Simplify则输出所求通解:⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧-++→])x 2[Sin ])1[c 4x (2]x 2[Cos ])2[c 81((e 81]x [y x 例2.4 (教材 例2.4) 求解微分方程x e x y +=''2, 并作出其积分曲线.输入g1=Table[Plot[E^x+x^3/3+c1+x*c2,{x,-5,5},DisplayFunction->Identity],{c1,-10,10,5},{c2,-5,5,5}];Show[g1,DisplayFunction->$DisplayFunction]; -4-224-40-20204060图2.3例2.5 (教材 例2.5) 求微分方程组⎪⎪⎩⎪⎪⎨⎧=--=++02y x dtdy e y x dt dx t 在初始条件0,100====t t y x 下的特解.输入122Clear[x,y,t];DSolve[{x' [t]+x[t]+2 y[t]==Exp[t], y'[t] -x[t]- y[t]==0,x[0]==1,y[0]==0},{x[t],y[t]},t]则输出所求特解:⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧+-→→])t [Sin ]t [Cos e (21]t [y ],t [Cos ]t [x t例2.6 验证c y y x =+--)3305(15152是微分方程2)(42-='y x x y 的通解. 输入命令<<Graphics`PlotField`<<Graphics`ImplicitPlot`sol=(-5x^3-30y+3y^5)/15==C;g1=ImplicitPlot[sol/.Table[{C->n},{n,-3,3}],{x,-3,3}];g2=PlotVectorField[{1,x^2/(y^4-2)},{x,-3,3},{y,-3,3},Frame->True,ScaleFunction->(1&),ScaleFactor->0.16,HeadLength->0.01,PlotPoints->{20,25}];g=Show[g2,g1,Axes->None,Frame->True];Show[GraphicsArray[{g1,g2,g}]];则分别输出积分曲线如图 2.4(a), 微分方程的方向场如图 2.4(b). 以及在同一坐标系中画出积分曲线和方向场的图形如下图2.4 (c).-3-2-1123-2-112-3-2-10123-3-2-10123-3-2-10123-3-2-10123图2.4从图 2.4(c)中可以看出微分方程的积分曲线与方向场的箭头方向吻合, 且当∞→x 时, 无论初始条件是什么, 所有的解都趋向于一条直线方程.例2.7 (教材 例2.6) 求解微分方程,)1(122/5+=+-x x y dx dy 并作出积分曲线. 输入<<Graphics`PlotField`DSolve[y' [x]-2y[x]/(x+1)==(x+1)^(5/2),y[x],x]则输出所给积分方程的解为 ⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧+++→]1[C )x 1()x 1(32]x [y 22/7123 下面在同一坐标系中作出这个微分方程的方向场和积分曲线(设),3,2,1,0,1,2,3---=C 输入t=Table[2(1+x)^(7/2)/3+(1+x)^2c,{c,-1,1}];g1=Plot[Evaluate[t],{x,-1,1},PlotRange->{{-1,1},{-2,2}},PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity];g2=PlotVectorField[{1,-2y/(x+1)+(x+1)^(5/2)},{x,-0.999,1},{y,-4,4},Frame->True,ScaleFunction->(1&), ScaleFactor->0.16,HeadLength->0.01,PlotPoints->{20,25},DisplayFunction->Identity];Show[g1,g2,Axes->None,Frame->True,DisplayFunction->$DisplayFunction];则输出积分曲线的图形（图2.5）.-0.75-0.5-0.2500.250.50.751-1.5-1-0.50.511.52图2.5例2.8 求解微分方程,2)21(22-+='-y x y xy 并作出其积分曲线.输入命令<<Graphics`PlotField`DSolve[1-2*x*y[x]*y' [x]==x^2+(y[x])^2-2,y[x],x]则得到微分方程的解为.)2(323C y x x y ++-+= 我们在33≤≤-C 时作出积分曲线, 输入命令t1=Table[(3+Sqrt[3])Sqrt[3+24x^2-4x^4-4*c*x]/(6*x),{c,-3,3}];t2=Table[(3-Sqrt[3])Sqrt[3+24x^2-4x^4-4*c*x]/(6*x),{c,-3,3}];gg1=Plot[Evaluate[t1],{x,-3,3},PlotRange->{{-3,3},{-3,3}},PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity];124gg2=Plot[Evaluate[t2],{x,-3,3},PlotRange->{{-3,3},{-3,3}},PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity];g1=ContourPlot[y-x^3/3-x*(-2+y^2),{x,-3,3},{y,-3,3},PlotRange->{-3,3},Contours->7,ContourShading->False,PlotPoints->50,DisplayFunction->Identity];g2=PlotVectorField[{1,(x^2+y^2-2)/(1-2*x*y)},{x,-3,3},{y,-3,3},Frame->True,ScaleFunction->(1&),ScaleFactor->0.16,HeadLength->0.01,PlotPoints->{20,25},DisplayFunction->Identity];Show[g1,g2,Axes->None,Frame->True,DisplayFunction->$DisplayFunction];Show[gg1,gg2,g2,Axes->None,Frame->True,DisplayFunction->$DisplayFunction];则输出微分方程的向量场与积分曲线, 并输出等值线的图2.6.-3-2-10123-2-10123-2-10123-2-1123图2.6用NDSolve 命令求微积分方程的近似解例2.9 (教材 例2.7) 求初值问题:1,0)1()1(2.1=='-++=x y y xy y xy 在区间[1.2,4]上的近似解并作图.输入fl=NDSolve[{(1+x*y[x])*y[x]+(1-x*y[x])*y'[x]==0,y[1.2]==1},y,{x,1.2,4}]则输出为数值近似解(插值函数)的形式:{{y->InterpolatingFunction[{{1.2,4.}},< >]}}用Plot 命令可以把它的图形画出来.不过还需要先使用强制求值命令Evalu-ate, 输入 Plot[Evaluate[y[x]/.fl],{x,1.2,4}]则输出近似解的图形（图2.7）.125 1.5 2.53 3.5410203040图2.7如果要求区间[1.2,4]内某一点的函数的近似值, 例如8.1=x y ,只要输入y[1.8]/.fl则输出所求结果{3.8341}例2.10 (教材 例2.8) 求范德波尔(Van der Pel)方程5.0,0,0)1(002-='==+'-+''==x x y y y y y y在区间[0,20]上的近似解.输入 Clear[x,y];NDSolve[{y''[x]+(y[x]^2-1)*y'[x]+y[x]==0,y[0]==0,y'[0]==-0.5},y,{x,0,20}];Plot[Evaluate[y[x]/.%],{x,0,20}]可以观察到近似解的图形（图2.8）.5101520-2-112图2.8126 ⎪⎩⎪⎨⎧==+-'1)1(01sin 2y x y x y x 的数值解, 并作出数值解的图形.输入命令<<Graphics`PlotField`sol=NDSolve[{x*y'[x]-x^2*y[x]*Sin[x]+1==0,y[1]==1},y[x],{x,1,4}];f[x_]=Evaluate[y[x]/.sol];g1=Plot[f[x],{x,1,4},PlotRange->All,DisplayFunction->Identity];g2=PlotVectorField[{1,(x^2*y*Sin[x]-1)/x},{x,1,4},{y,-2,9},Frame->True,ScaleFunction->(1&),ScaleFactor->0.16,HeadLength->0.01,PlotPoints->{20,25},DisplayFunction->Identity];g=Show[g1,g2,Axes->None,Frame->True];Show[GraphicsArray[{g1,g}],DisplayFunction->$DisplayFunction];则输出所给微分方程的数值解及数值解的图2.9.1.522.533.544681 1.52 2.53 3.54-22468例2.11 (教材 例2.9) 求出初值问题⎪⎩⎪⎨⎧='==+'+''0)0(,1)0(cos sin 22y y xy x y y的数值解, 并作出数值解的图形.输入NDSolve[{y''[x]+Sin[x]^2*y'[x]+y[x]==Cos[x]^2,y[0]==1,y'[0]==0},y[x],{x,0,10}]127 Plot[Evaluate[y[x]/.%],{x,0,10}];则输出所求微分方程的数值解及数值解的图形（图2.10）.2468100.20.40.60.8图2.10例2.12 (教材 例2.10) 洛伦兹(Lorenz)方程组是由三个一阶微分方程组成的方程组.这三个方程看似简单, 也没有包含复杂的函数, 但它的解却很有趣和耐人寻味. 试求解洛伦兹方程组,0)0(,4)0(,12)0()(4)()()()()(45)()()()(16)(16)(⎪⎪⎩⎪⎪⎨⎧===-='-+-='-='z y x t z t y t x t z t y t x t z t x t y t x t y t x 并画出解曲线的图形.输入Clear[eq,x,y,z]eq=Sequence[x'[t]==16*y[t]-16*x[t],y'[t]==-x[t]*z[t]-y[t]+45x[t],z'[t]==x[t]*y[t]-4z[t]];sol1=NDSolve[{eq,x[0]==12,y[0]==4,z[0]==0},{x[t],y[t],z[t]},{t,0,16},MaxSteps->10000];g1=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.sol1],{t,0,16},PlotPoints->14400,Boxed->False,Axes->None];则输出所求数值解的图形（图2.11(a)）. 从图中可以看出洛伦兹微分方程组具有一个奇异吸引子, 这个吸引子紧紧地把解的图形“吸”在一起. 有趣的是, 无论把解的曲线画得多长, 这些曲线也不相交.128图2.11改变初值为,10)0(,10)0(,6)0(=-==z y x 输入sol2=NDSolve[{eq,x[0]==6,y[0]==-10,z[0]==10}, {x[t],y[t],z[t]},{t,0,24},MaxSteps->10000];g2=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.sol2],{t,0,24},PlotPoints->14400,Boxed->False,Axes->None];Show[GraphicsArray[{g1,g2}]];则输出所求数值解的图形（图2.11(b)）. 从图中可以看出奇异吸引子又出现了, 它把解“吸”在某个区域内, 使得所有的解好象是有规则地依某种模式缠绕.实验习题1. 求下列微分方程的通解:(1) ;0136=+'+''y y y(2) ();024=+''+y y y(3) ;2sin 52x e y y y x =+'-''(4) .)1(963x e x y y y +=+'-''2. 求下列微分方程的特解:(1) ;15,0,029400='==+'+''==x x y y y y y(2) .1,1,02sin ='==++''==ππx x y yx y y 3. 求微分方程0cos 2)1(2=-+'-x xy y x 在初始条件10==x y 下的特解.分别求精确解和数值解)10(≤≤x 并作图.4. 求微分方程组⎪⎪⎩⎪⎪⎨⎧=--=++t t e y x dt dy e y x dt dx 235的通解.129 5. 求微分方程组⎪⎪⎩⎪⎨⎧==+-==-+==4,081,0300t t y y x dt dyxy x dt dx 的特解. 6. 求欧拉方程组324x y y x y x =-'+''的通解.7. 求方程5,0,011='==+'+''==x x y y y y x y 在区间[0,4]上的近似解.。