05Mathematica项目四 无穷级数与微分方程

实验1 无穷级数（基础实验）

实验目的

观察无穷级数部分和的变化趋势,进一步理解级数的审敛法以及幂级数部分和对函数的 逼近. 掌握用Mathematica求无穷级数的和, 求幂级数的收敛域, 展开函数为幂级数以及展 开周期函数为傅里叶级数的方法.

基本命令

1. 求无穷和的命令Sum

该命令可用来求无穷和. 例如,输入 Sum[1/n^2,{n,l,Infinity}]

则输出无穷级数的和为2/6. 命令Sum与数学中的求和号相当. 2. 将函数展开为幂级数的命令Series 该命令的基本格式为

Series[f[x],{x,x0,n}]

它将f(x)展开成关于xx0的幂级数. 幂级数的最高次幂为(xx0)n,余项用(xx0)n1表 示. 例如,输入

Series[y[x],{x,0,5}] 则输出带皮亚诺余项的麦克劳林级数

1114156y0x4y0xox y0y0xy0x2y30x326241203. 去掉余项的命令Normal

在将f(x)展开成幂级数后, 有时为了近似计算或作图, 需要把余项去掉. 只要使用 Normal命令. 例如,输入

Series[Exp[x],{x,0,6}] Normal[%] 则输出

x2x3x4x5x6O[x]7 1x2!3!4!5!6!x2x3x4x5x6 1x 2!3!456! 4. 强制求值的命令Evaluate

如果函数是用Normal命令定义的, 则当对它进行作图或数值计算时, 可能会出现问题. 例如,输入

fx=Normal[Series[Exp[x],{x,0,3}]] Plot[fx,{x,-3,3}]

则只能输出去掉余项后的展开式

x2x3 1x 26而得不到函数的图形. 这时要使用强制求值命令Evaluate, 改成输入 Plot[Evaluate[fx],{x,-3,3}] 则输出上述函数的图形.

5. 作散点图的命令ListPlot

119

ListPlot [ ]为平面内作散点图的命令, 其对象是数集,例如,输入

ListPlot[Table[j^2,{j,16}],PlotStyle->PointSize[0,012]]

则输出坐标为{1,12},{2,22},{3,32},,{16,162}的散点图(图1.1).

250200150100502.557.51012.515 图1.1

6. 符号“/；”用于定义某种规则,“/；”后面是条件. 例如,输入

Clear[g,gf];

g[x_]:=x/;0<=x<1 g[x_]:=-x/;-1<=x<0 g[x_]:=g[x –2]/;x>=1

则得到分段的周期函数

1x0x,g(x)x,0x1

g(x2),x1再输入

gf=Plot[g[x],{x,-1,6}] 则输出函数g(x)的图形1.2.

10.8 0.60.40.2-1123456 图1.2

注:用Which命令也可以定义分段函数, 从这个例子中看到用“…(表达式)/; …(条件)”来 定义周期性分段函数更方便些. 用Plot命令可以作出分段函数的图形, 但用Mathematica命 令求分段函数的导数或积分时往往会有问题. 用Which定义的分段函数可以求导但不能积 120

分. Mathematica内部函数中有一些也是分段函数. 如:Mod[x,1],Abs[x],Floor[x]和UnitStep[x]. 其中只有单位阶跃函数UnitStep[x]可以用Mathematica命令来求导和求定积分. 因此在求分 段函数的傅里叶系数时, 对分段函数的积分往往要分区来积. 在被积函数可以用单位阶跃函 数UnitStep的四则运算和复合运算表达时, 计算傅里叶系数就比较方便了.

实验举例 数项级数

例1.1 (教材 例1.1)

(1) 观察级数12的部分和序列的变化趋势.

n1n(2) 观察级数1n的部分和序列的变化趋势.

n1输入

s[n_]=Sum[1/k^2,{k,n}];data=Table[s[n],{n,100}]; ListPlot[data];

N[Sum[1/k^2,{k,Infinity}]] N[Sum[1/k^2,{k,Infinity}],40]

则输出(1)中级数部分和的变化趋势图1.3. 1.62204060801001.581.561.54图1.3 级数的近似值为1.64493.

输入

s[n_]=Sum[1/k,{k,n}];data=Table[s[n],{n,50}]; ListPlot[data,PlotStyle->PointSize[0.02]];

则输出(2)中级数部分和的的变化趋势图1.4. 4.543.532.521.51020304050图1.4

121

例1.2 (教材 例1.2) 画出级数

(1)n1n11的部分和分布图. n输入命令

Clear[sn,g];sn=0;n=1;g={};m=3;

While[1/n>10^-m,sn=sn+(-1)^(n-1)/n;

g=Append[g,Graphics[{RGBColor[Abs[Sin[n]],0,1/n],

Line[{{sn,0},{sn,1}}]}]];n++];

Show[g,PlotRange->{-0.2,1.3},Axes->True];

则输出所给级数部分和的图形（图1.5）,从图中可观察到它收敛于0.693附近的一个数.

1.210.80.60.40.20.688-0.20.6920.6940.6960.698图1.5 例1.3 求

1的值. 24n8n3n1输入

Sum[x^(3k),{k,1,Infinity}] 得到和函数

x3 

1x3例1.4 (教材 例1.3) 设an 输入

Clear[a];

a[n_]=10^n/(n!);

vals=Table[a[n],{n,1,25}];

ListPlot[vals,PlotStyle->PointSize[0.012]]

则输出an的散点图（1.6）,从图中可观察an的变化趋势. 输入 Sum[a[n],{n,l,Infinity}] 则输出所求级数的和.

10n, 求n!an1n.

122

2500200015001000500510152025 图1.6

求幂级数的收敛域 例1.5 (教材 例1.4) 求 输入

Clear[a];

a[n_]=4^(2n)*(x-3)^n/(n+1); stepone=a[n+1]/a[n]//Simplify

则输出

16(1n)(3x)

2n n042n(x3)n的收敛域与和函数.

n1再输入

steptwo=Limit[stepone,n->Infinity] 则输出

16(3x)

这里对a[n+1]和a[n]都没有加绝对值. 因此上式的绝对值小于1时, 幂级数收敛; 大于1 时发散. 为了求出收敛区间的端点, 输入

ydd=Solve[steptwo==1,x] zdd=Solve[steptwo==-1,x]

则输出

4947 x与x

161647494749x时,级数收敛,当x或x时,级数发散.

16161616为了判断端点的敛散性, 输入 Simplify[a[n]/.x->(49/16)] 则输出右端点处幂级数的一般项为

1 n149因此,在端点x处,级数发散. 再输入

16 Simplify[a[n]/.x->(47/16)]

由此可知,当

123

则输出左端点处幂级数的一般项为

(1)n n147因此,在端点x处, 级数收敛.

16也可以在收敛域内求得这个级数的和函数. 输入 Sum[4^(2n)*(x-3)^n/(n+1),{n,0,Infinity}] 则输出

Log[116(3x)] 

16(3x)

函数的幂级数展开

例1.6 (教材 例1.5) 求cosx的6阶麦克劳林展开式. 输入

Series[Cos[x],{x,0,6}] 则输出

x2x4x6o[x]7 1224720注:这是带皮亚诺余项的麦克劳林展开式.

例1.6 (教材 例1.6) 求lnx在x1处的6阶泰勒展开式.

输入

Series[Log[x],{x,1,6}] 则输出

(x1)2(x1)3(x1)4(x1)5(x1)6o[x]7. 23456例1.7 (教材 例1.7) 求arctanx的5阶泰勒展开式. 输入

serl=Series[ArcTan[x],{x,0,5}]; Poly=Normal[serl]

则输出arctanx的近似多项式

(x1)x3x5 x 35 通过作图把arctanx和它的近似多项式进行比较. 输入

Plot[Evaluate[{ArcTan[x],Poly}],{x,-3/2,3/2},

PlotStyle->{Dashing[{0.01}],GrayLevel[0]},AspectRatio->l]

则输出所作图形（图1.7）, 图中虚线为函数arctanx,实线为它的近似多项式.

124

1.510.5-1.5-1-0.5-0.50.511.5-1-1.5图1.7

x12x12 例1.9 求e在x1处的8阶泰勒展开, 并通过作图比较函数和它的近似多项式. 输入

Clear[f];

f[x_]=Exp[-(x-1)^2*(x+1)^2];

poly2=Normal[Series[f[x],{x,1,8}]]

Plot[Evaluate[{f[x],poly2}],{x,-1.75,1.75},PlotRange->

{-2,3/2},PlotStyle->{Dashing[{0.01}], GrayLevel[0]}]

则得到近似多项式和它们的图1.8.

141x41x71x161x

2345

41x6281x71731x8 361.510.5-1.5-1-0.5-0.5-1-1.5-20.511.5 图1.8

125

例1.10 求函数sinx在x0处的3,5,7,,91阶泰勒展开, 通过作图比较函数和它的近似多项式, 并形成动画进一步观察. 因为

x2k11o(x2n2) sinx2k1!k0nk所以输入

Do[Plot[{Sum[(-1)^j*x^(2j+1)/(2j+1)!,{j,0,k}],

Sin[x]},{x,-40,40},

PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]}],{k,1,45}]

则输出为sinx的3阶和91阶泰勒展开的图形. 选中其中一幅图形,双击后形成动画. 图1.9是最后一幅图.

42-40-202040-2-4图1.9

例1.11 利用幂级数展开式计算5240(精确到1010). 因为

5 124052433314.

315根据(1x)m在x0处的展开式有

5141149111240314283 125352!353!3k15i111n11i1Sn314k34k53k25k!故前n(n2)项部分和为

 输入命令

s[n_]=3(1-1/(5*3^4)-Sum[Product[5i-1,{i,1,k-1}]/(5^k k!3^(4k)),{k,2,n-1}]); r[n_]=Product[5i-1,{i,1,n-1}]/5^n/n!3^(4n-5)/80; 126

delta=10^(-10);n0=100;

Do[Print[\"n=\

If[r[n]则输出结果为

傅里叶级数

例1.12 (教材 例1.8) 设g(x)是以2为周期的周期函数,它在,的表达式是

52402.992555739.

1,x0 g(x)

1,0x将g(x)展开成傅里叶级数. 输入

Clear[g];

g[x_]:=-1/;-Pi<=x<0 g[x_]:=1/;0<=xPlot[g[x],{x,-Pi,5 Pi},PlotStyle->{RGBColor[0,1,0]}];

则输出g(x)的图形 (图1.10).

10.5-2.52.557.51012.515-0.5-1图1.10

因为g(x)是奇函数, 所以它的傅里叶展开式中只含正弦项. 输入

b2[n_]:=b2[n]=2 Integrate[1*Sin[n*x],{x,0,Pi}]/Pi; fourier2[n_,x_]:=Sum[b2[k]*Sin[k*x],{k,1,n}];

tu[n_]:=Plot[{g[x],Evaluate[fourier2[n,x]]}, {x,-Pi,5 Pi},

PlotStyle->{RGBColor[0,1,0],RGBColor[1,0.3,0.5]},DisplayFunction->Identity];

(*tu[n]是以n为参数的作图命令*)

tu2=Table[tu[n],{n,1,30,5}];

(*tu2是用Table命令作出的6个图形的集合*)

toshow=Partition[tu2,2];

(*Partition是对集合tu2作分割, 2为分割的参数*)

Show[GraphicsArray[toshow]]

(*GraphicsArray是把图形排列的命令*)

则输出6个排列着的图形（图1.11）,每两个图形排成一行. 可以看到n越大, g(x)的傅里叶级数的前

n项和与g(x)越接近.

127

10.5-2.5-0.5-110.5-2.5-0.5-110.5-2.5-0.5-12.557.51012.5152.557.51012.5152.557.51012.51510.5-2.5-0.5-110.5-2.5-0.5-110.5-2.5-0.5-12.557.51012.5152.557.51012.5152.557.51012.515 图1.11

实验习题

1.求下列级数的和:

(1)



k1



k

; (2)2k

(2k1)k112; (3)

(2k)k112; (4)

k1(1)k1. k 2. 求幂级数

n0(x1)2n1(5)n的收敛域与和函数.

3. 求函数(1x)ln(1x)的6阶麦克劳林多项式. 4. 求arcsinx的6阶麦克劳林多项式.

x5. 设f(x)2,求f(x)的5阶和10阶麦克劳林多项式,把两个近似多项式和函数的

x1图形作在一个坐标系内.

116. 设f(x)在一个周期内的表达式为f(x)1x2x, 将它展开为傅里叶级

22数(取6项), 并作图.

0x11,7. 设f(x)在一个周期内的表达式为f(x), 将它展开为傅里叶级数

2x,1x2(取8项), 并作图. 8. 求级数 128

sink的和的近似值. kk1

实验2 微分方程（基础实验）

实验目的 理解常微分方程解的概念以及积分曲线和方向场的概念,掌握利用 Mathematica求微分方程及方程组解的常用命令和方法.

基本命令

1. 求微分方程的解的命令DSolve

对于可以用积分方法求解的微分方程和微分方程组,可用Dsolve命令来求其通解或特解. 例如,求方程y3y2y0的通解, 输入

DSolve[y ''[x]+3y '[x]+2y[x]==0,y[x],x]

则输出含有两个任意常数C[1]和C[2]的通解:

y[x]e2xC[1]exC[2]

注:在上述命令中,一阶导数符号 ' 是通过键盘上的单引号 ' 输入的,二阶导数符号 '' 要 输入两个单引号,而不能输入一个双引号.

又如,求解微分方程的初值问题:

y4y3y0,yx06,yx010,

输入

Dsolve[{y''[x]+4 y'[x]+3y[x]==0,y[0]==6, y'[0]==10},y[x],x]

(*大括号把方程和初始条件放在一起*) 则输出

y[x]e3x(814e2x

2. 求微分方程的数值解的命令NDSolve

对于不可以用积分方法求解的微分方程初值问题,可以用NDSolve命令来求其特解.例如 要求方程

yy2x3,yx00.5 的近似解(0x1.5), 输入

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命令得到的输出是解yy(x)的近似值. 首先在区间[0,1.5]内插入一系 列点x1,x2,,xn, 计算出在这些点上函数的近似值y1,y2,,yn, 再通过插值方法得到

yy(x)在区间上的近似解.

3. 一阶微分方程的方向场

一般地,我们可把一阶微分方程写为

yf(x,y)

的形式,其中f(x,y)是已知函数. 上述微分方程表明:未知函数y在点x处的斜率等于函数

f在点(x,y)处的函数值. 因此,可在Oxy平面上的每一点, 作出过该点的以f(x,y)为斜率 的一条很短的直线(即是未知函数y的切线). 这样得到的一个图形就是微分方程yf(x,y) 的方向场. 为了便于观察, 实际上只要在Oxy平面上取适当多的点,作出在这些点的函数的 切线. 顺着斜率的走向画出符合初始条件的解,就可以得到方程yf(x,y)的近似的积分曲 线.

129

例如, 画出输入

dy1y2,y(0)0的方向场. dx<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）, 从图中可以观察到, 当初始条件为y01/2时, 这个微分方程的解介于1和1之间, 且当x趋向于或时, y(x)分别趋向于1与1.

221100-1-1-2-3-2-10123-2-2-10123 -3 图2.1 图2.2

下面求解这个微分方程, 并在同一坐标系中画出方程的解与方向场的图解. 输入

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];

1e2x则输出微分方程的解y(x),以及解曲线与方向场的图形（图2.2）. 从图中可以看

1e2x到, 微分方程的解与方向场的箭头方向相吻合.

实验内容

用Dsolve命令求解微分方程

例2.1 (教材 例2.1) 求微分方程 y2xyxex的通解. 输入

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] 则输出微分方程的通解:

221 y[x]exx2exC[1]

2其中C[1]是任意常数.

例2.2 (教材 例2.2) 求微分方程xyyex0在初始条件y 输入

Clear[x,y];

DSolve[{x*y ' [x]+y[x]-Exp[x]==0,y[1]==2 E},y[x],x]

则输出所求特解: 130

x122e下的特解.

eex y[x]x  例2.3 (教材 例2.3) 求微分方程y2y5yexcos2x的通解.

输入

DSolve[y ''[x]-2y '[x]+5y[x]==Exp[x]*Cos[2 x],y[x],x]//Simplify

则输出所求通解:

1 y[x]ex((18c[2])Cos[2x]2(x4c[1])Sin[2x])

8例2.4 (教材 例2.4) 求解微分方程y2xex, 并作出其积分曲线.

输入

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];

则输出积分曲线的图形（图2.3）. 604020-4-2-20-4024 图2.3

dxx2yetdt例2.5 (教材 例2.5) 求微分方程组在初始条件xt01,yt00下的特

dyxy0dt解.

输入

Clear[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]

则输出所求特解:

1 x[t]Cos[t],y[t](etCos[t]Sin[t])

2

1x225(5x30y3y)c是微分方程y(x)4例2.6 验证的通解. 15y2输入命令

131

<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).

21-3-2-1-1-21233210-1-2-3-3-2-101233210-1-2-3-3-2-10123(a) (b) (c) 图2.4

从图2.4(c)中可以看出微分方程的积分曲线与方向场的箭头方向吻合, 且当x时, 无论初始条件是什么, 所有的解都趋向于一条直线方程.

例2.7 (教材 例2.6) 求解微分方程输入

<DSolve[y' [x]-2y[x]/(x+1)==(x+1)^(5/2),y[x],x]

则输出所给积分方程的解为

dy2y(x1)5/2,并作出积分曲线. dxx12 y[x](1x)7/2(1x)2C[1]

3下面在同一坐标系中作出这个微分方程的方向场和积分曲线(设C3,2,1,0,1,2,3), 输入

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];

132

则输出积分曲线的图形（图2.5）.

21.510.50-0.5-1-1.5-0.75-0.5-0.2500.250.50.751 图2.5

例2.8 求解微分方程(12xy)yx2y22,并作出其积分曲线. 输入命令

<DSolve[1-2*x*y[x]*y' [x]==x^2+(y[x])^2-2,y[x],x]

则得到微分方程的解为

x3yx(2y2)C.

3 我们在3C3时作出积分曲线, 输入命令

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];

gg2=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,Cont

ourShading->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,

133

DisplayFunction->$DisplayFunction];

则输出微分方程的向量场与积分曲线, 并输出等值线的图2.6.

3213021-10-1-2-2-3-2-10123

-2-10123 图2.6

用NDSolve命令求微积分方程的近似解

例2.9 (教材 例2.7) 求初值问题:(1xy)y(1xy)y0,yx1.21在区间[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）.

403020101.52.533.54 图2.7

如果要求区间[1.2,4]内某一点的函数的近似值, 例如y y[1.8]/.fl 则输出所求结果

{3.8341}

134

x1.8 ,只要输入

例2.10 (教材 例2.8) 求范德波尔(Van der Pel)方程

y(y21)yy0,yx00,yx00.5

在区间[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）.

215101520-1-2 图2.8

xyx2ysinx10y(1)1 的数值解, 并作出数值解的图形.

输入命令

<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];

135

g=Show[g1,g2,Axes->None,Frame->True]; Show[GraphicsArray[{g1,g}],DisplayFunction->

$DisplayFunction];

则输出所给微分方程的数值解及数值解的图2.9.

8641.5

例2.11 (教材 例2.9) 求出初值问题

22.533.5486420-211.522.533.54 图2.9 22yysinxycosx y(0)1,y(0)0的数值解, 并作出数值解的图形.

输入

NDSolve[{y''[x]+Sin[x]^2*y'[x]+y[x]==Cos[x]^2,

y[0]==1,y'[0]==0},y[x],{x,0,10}]

Plot[Evaluate[y[x]/.%],{x,0,10}];

则输出所求微分方程的数值解及数值解的图形（图2.10）.

20.8468100.60.40.2 图2.10

例2.12 (教材 例2.10) 洛伦兹(Lorenz)方程组是由三个一阶微分方程组成的方程组.这三个方程看似简单, 也没有包含复杂的函数, 但它的解却很有趣和耐人寻味. 试求解洛伦兹方程组 136

x(t)16y(t)16x(t)y(t)x(t)z(t)45x(t)y(t), z(t)x(t)y(t)4z(t)x(0)12,y(0)4,z(0)0并画出解曲线的图形.

输入

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)）. 从图中可以看出洛伦兹微分方程组具有一个奇异吸引子, 这个吸引子紧紧地把解的图形“吸”在一起. 有趣的是, 无论把解的曲线画得多长, 这些曲线也不相交.

(a) (b) 图2.11

改变初值为x(0)6,y(0)10,z(0)10,输入

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)）. 从图中可以看出奇异吸引子又出现了, 它把解“吸”在某个区域内, 使得所有的解好象是有规则地依某种模式缠绕.

137

实验习题

1. 求下列微分方程的通解:

(1) y6y13y0;

(2) y42yy0; (3) y2y5yexsin2x;

(4) y6y9y(x1)e3x. 2. 求下列微分方程的特解:

(1) y4y29y0,yx00,yx015; (2) yysin2x0,y和数值解(0x1)并作图.

dxtdt5xye 4. 求微分方程组的通解.

dy2tx3yedtdx3xy0,xt01dt 5. 求微分方程组的特解.

dy8xy0,yt04dtx1,yx1.

x03. 求微分方程(x21)y2xycosx0在初始条件y1下的特解.分别求精确解

6. 求欧拉方程组x2yxy4yx3的通解.

7. 求方程yxyy0,yx10,yx15在区间[0,4]上的近似解.

实验3 抛射体的运动（续）（综合实验）

实验目的 通过微分方程建模和Mathematica软件,在项目一实验5的基础上,进一步研 究在考虑空气阻力的情况下抛射体的运动.

问题 根据侦察,发现离我军大炮阵地水平距离10km的前方有一敌军的坦克群正以每小 时50km向我军阵地驶来,现欲发射炮弹摧毁敌军坦克群. 为在最短时间内有效摧毁敌军坦 克,要求每门大炮都能进行精射击,这样问题就可简化为单门大炮对移动坦克的精确射击 问题. 假设炮弹发射速度可控制在0.2km/s至0.5km/s之间,问应选择怎样的炮弹发射速度和 怎样的发射角度可以最有效摧毁敌军坦克.

说明 本节我们研究受到重力和空气阻力约束的抛射体的射程. 用(x(t),y(t))记抛射体 的位置, 其中x轴是运动的水平方向, y轴是垂直方向. 通过在y0的约束下最大化x, 可以 计算出使抛射体的射程最大的发射角. 假设t0时抛射体(炮弹)在原点(0,0)以与水平线夹角 为,初始速度为v0发射出去. 它受到的空气阻力为

dxdy Frkvk,. (3.1)

dtdt重力为 138

Fg(0,mg). (3.2) 在推导x(t)和y(t)所满足的微分方程之前, 补充一点说明:虽然我们将位置变量

x(t),y(t)仅写作t的函数,但实际上位置变量还依赖于几个其它的变量或参数. 特别是,x和y 也依赖于发射角、阻力系数k、质量m及重力加速度g等.

为了推导x和y的方程, 按照牛顿定律Fma,并结合重力的公式(3.2)和空气阻力的公 式(3.1), 得到微分方程

mx(t)kx(t)0 (3.3) my(t)ky(t)mg0 (3.4)

根据前面所述假设知, x(t),y(t)满足下列初始条件

x(0)0,y(0)0,x(0)v0cos,y(0)v0sin. (3.5)

先求解x(t),由方程(3.3),令vx,可将其化为一阶微分方程

mvkv0.

易求出其通解 v(t)由v(0)x(0)v0cos, 得Cv0cos,所以

ktmCe.

v(t)v0从vx,通过积分得到x, 即

ktmcose.

tmx(t)v0cosemD.

kkm由x(0)0, 得Dv0cos, 所以

ktmx(t)v0cos(1em) (3.6)

k类似地,可从方程(3.4)解出y. 令vy, 方程化为一阶微分方程, 两端除以m,得

kvvg.

mkktm再在上述方程两端乘以积分因子e.ktemv得

kktkmtevgem, mttd即 (vem)gem,

dtkk两端积分得

ktvemgmmteC.

kkktgm所以 vCem.

k利用初始条件y(0)v(0)v0sin确定其中的常数C后, 积分v得到y,再次利用初始条

件y(0)0确定任意常数后,则得到

139

tgmmmmtmv(1em)sin. (3.7) yte0kkkk下面我们利用公式(3.6)与(3.7)来描绘炮弹运行的典型图形.

kk假定炮弹发射的初速度为0.25km/s, 发射角为55, 输入

Clear[a,t,x,y,g,m,k]

x[v_,a_,t_]:=(m/k)*v*Cos[a Pi/180]*(1-Exp[-(k/m)*t]) y[v_,a_,t_]:=(g*m/k)((m/k)-t-(m/k)*Exp[-(k/m)*t])+

(m/k)*v*Sin[a Pi/180]*(1-Exp[-(k/m)*t])

g=9.8;m=5.0;k=0.01;

炮弹飞行的时间由炮弹落地时的条件y0所确定. 输入 FindRoot[y[350,55,t]==0,{t,50}]

则输出炮弹飞行的时间 {t->57.4124}

当发射角65时, 输入

x[350,55, 57.4124]//N 则输出炮弹的最大射程为 10888.5

现在我们可以画出炮弹运行的典型轨迹了. 输入 ParametricPlot[{x[350,55,t],y[350,55,t]},

{t,0,57.4124},PlotRange->{0,11000},AxesLabel->{x,y}]

则输出图3.1. y100008000600040002000x200040006000800010000

图3.1

实验报告

在上述假设下,进一步研究下列问题:

(1) 选择一个初始速度和发射角,利用Mathematica画出炮弹运行的典型轨迹.

(2) 假定坦克在大炮前方10km处静止不动,炮弹发射的初速度为0.32km/s,应选择什么样 的发射角才能击中坦克?画出炮弹运行的几个轨迹图,通过实验数据和图形来说明你的结论 的合理性.

(3) 假定坦克在大炮前方10km处静止不动,探索降低或调高炮弹发射的初速度的情况 下,应如何选择炮弹的发射角?从上述讨论中总结出最合理有效的发射速度和发射角.

(4) 在上题结论的基础上,继续探索,假定坦克在大炮前方10km处以每小时50km向大炮 方向前进,此时应如何制定迅速摧毁敌军坦克的方案? 140

注:在研究过程中,还要包括适当改变阻力系数k与炮弹的质量m所带来的变化.

实验4 蹦极跳运动（综合实验）

实验目的 利用Mathematica软件,通过微分方程建模,研究蹦极跳运动.

问题 在不考虑空气阻力和考虑空气阻力等多种情况下,研究蹦极跳运动中,蹦极者与蹦 极绳设计之间的各种关系.

说明 蹦极绳相当于一根粗橡皮筋或有弹性的绳子. 当受到张力使之超过其自然长度,绳 子会产生一个线性回复力, 即绳子会产生一个力使它恢复到自然长度, 而这个力的大小与它 被拉伸的长度成正比. 在一次完美的蹦极跳过程中, 蹦极者爬上一座高桥或高的建筑物, 把 绳的一头系在自己身上, 另一头系在一个固定物体如桥栏杆上, 当他跳离桥时, 激动人心的 时刻就到来了. 这里要分析的是蹦极者从跳出那一瞬间起他的运动规律.

首先要建立坐标系. 假设蹦极者的运动轨迹是垂直的, 因此我们只要用一个坐标来确 定他在时刻t的位置. 设y是垂直坐标轴, 单位为英尺, 正向朝下, 选择y0为桥平面, 时间 t的单位为秒, 蹦极者跳出的瞬间为t0, 则y(t)表示t时刻蹦极者的位置. 下面我们要求出

y(t)的表达式.

由牛顿第二定律, 物体的质量乘以加速度等于物体所受的力. 我们假设蹦极者所受的力 只有重力、空气阻力和蹦极绳产生的回复力. 当然, 直到蹦极者降落的距离大于蹦极绳的自 然长度时, 蹦极绳才会产生回复力. 为简单起见, 假设空气阻力的大小与速度成正比, 比例 系数为1, 蹦极绳回复力的比例系数为0.4. 这些假设是合理的, 所得到的数学结果与研究所 做的蹦极实验非常吻合. 重力加速度g32ft/s2.

现在我们来考虑一次具体的蹦极跳. 假设绳的自然长度为L200ft, 蹦极者的体重为 ①②

160lb,则他的质量为m160/325斯. 在他到达绳的自然长度(即yL200)前, 蹦 极者的坠落满足下列初值问题:

dy1gv, v(0)0. dtm利用Mathematica求解上述问题. 输入

g=32; m=5; L=200;

{{v1[t_],y1[t_]}}={v[t],y[t]}/.

DSolve[{v'[t]==-g-v[t]/m,y'[t]==v[t],v[0]==0,y[0]==0},{v,y},t]

则输出

{{160et/5(1et/5),160et/5(55et/5et/5t)}}

蹦极者坠落L英尺所用的时间为

t1=t/.FindRoot[y1[t]==-L,{t,2}] 4.00609

现在我们需要找到当蹦极绳产生回复力后的运动初始条件. 当tt1时, 蹦极者的坠落 满足方程

dv10.4gv(Ly) dtmm初始条件为y(t1)L,v(t1)v1(t1).解初值问题:

{{v2[t_],y2[t_]}}={v[t],y[t]}/.

DSolve[{v'[t]==-g-v[t]/m-0.4*(L+y[t])/m,

y'[t]==v[t],v[t1]==v1[t1],y[t1]==-L},{v,y},t]

141

则输出下列结果

i)t{{(0.1469210.262116i)e(0.10.264575((299.96173.3019i)(4.0167310141119.99i)(0.0.52915i)t2.11982i)(0.0.52915i)te0.400609(2.6855710141.074231013i)e(0.400609_1.05991i)(0.0.52915i)t(617.364233.342i)e(0.801218),i)t (0.1469210.262116i)e(0.10.264575((132.5281083.65i)(3704.011399.98i)(0.0.52915i)t1.05991i)(0.0.52915i)te0.400609(0.2333.42i)e(0.801218_(0.10.264575i)t(4200.1587.45i)e0.400609_2.11982i)(0.10.264575i)t(4200.1587.45i)e(0.400609)}}

这个解是用复指数函数来表示的.

现在蹦极者的位置由命令

bungeey[t_]=If[t给出, 输入命令

Plot[bungeey[t],{t,0,40},PlotRange->All]

则输出位置-时间图形(图4.1)

10-200203040-400-600

图4.1

从上图可以看出, 蹦极者在大约13s内由桥面坠落770ft, 然后弹回到桥面下550ft, 上下 振动几次, 最终降落到桥面下大约600ft处.

实验报告

1.在上述问题中(L200,w160),求出需要多长时间蹦极者才能到达他运动轨迹上的 最低点, 他能下降到桥面下多少英尺?

2.用图描述一个体重为195lb, 用200ft长绳子的蹦极者的坠落. 在绳子对他产生力之前, 他能做多长时间的“自由”降落?

3.假设你有一根300ft长的蹦极索, 在一组坐标轴上画出你所在实验组的全体成员的运 动轨迹草图.

4.一个55岁, 体重185lb的蹦极者, 用一根250ft长的蹦极索. 在降落过程中, 他达到的 最大速度是多少? 当他最终停止运动时, 他被挂在桥面下多少英尺?

5.用不同的空气阻力系数和蹦极索常数做实验, 确定一组合理的参数, 使得在这组参数 142

下, 一个160lb的蹦极者可以回弹到蹦极索的自然长度以上.

6.科罗拉多的皇家乔治桥(它跨越皇家乔治峡谷)距谷底1053ft, 一个175lb的蹦极者希望 能正好碰到谷底, 则他应使用多长的绳子?

7.假如上题中的蹦极者体重增加10lb, 再用同样长的绳子从皇家乔治桥上跳下, 则当他 撞到乔治峡谷谷底时, 他的坠落速度是多少?

143