项目四无穷级数与微分方程.doc

项目四 无穷级数与微分方程 实验1 无穷级数（基础实验） 实验目的观察无穷级数部分和的变化趋势,进一步理解级数的审敛法以及幂级数部分和对函数的逼近. 掌握用Mathematica求无穷级数的和, 求幂级数的收敛域, 展开函数为幂级数以及展开周期函数为傅里叶级数的方法. 数项级数例1.1 (教材 例1.1)(1) 观察级数的部分和序列的变化趋势.(2) 观察级数的部分和序列的变化趋势.输入sn_=Sum1/k2,k,n;data=Tablesn,n,100;ListPlotdata;NSum1/k2,k,InfinityNSum1/k2,k,Infinity,40则输出(1)中级数部分和的变化趋势图1.3.图1.3级数的近似值为1.64493.输入sn_=Sum1/k,k,n;data=Tablesn,n,50;ListPlotdata,PlotStyle-PointSize0.02;则输出(2)中级数部分和的的变化趋势图1.4.图1.4例1.2 (教材 例1.2) 画出级数的部分和分布图.输入命令Clearsn,g;sn=0;n=1;g=;m=3;While1/n10-m,sn=sn+(-1)(n-1)/n;g=Appendg,GraphicsRGBColorAbsSinn,0,1/n,Linesn,0,sn,1;n+;Showg,PlotRange-0.2,1.3,Axes-True;则输出所给级数部分和的图形（图1.5）,从图中可观察到它收敛于0.693附近的一个数.图1.5例1.3 求的值.输入 Sumx(3k),k,1,Infinity得到和函数 例1.4 (教材 例1.3) 设 求. 输入Cleara;an_=10n/(n!);vals=Tablean,n,1,25;ListPlotvals,PlotStyle-PointSize0.012则输出的散点图（1.6）,从图中可观察的变化趋势. 输入 Suman,n,l,Infinity则输出所求级数的和.图1.6求幂级数的收敛域 例1.5 (教材 例1.4) 求的收敛域与和函数. 输入Cleara;an_=4(2n)*(x-3)n/(n+1);stepone=an+1/an/Simplify则输出 再输入 steptwo=Limitstepone,n-Infinity则输出 这里对an+1和an都没有加绝对值. 因此上式的绝对值小于1时, 幂级数收敛; 大于1时发散. 为了求出收敛区间的端点, 输入ydd=Solvesteptwo=1,xzdd=Solvesteptwo=-1,x则输出 由此可知,当时,级数收敛,当或时,级数发散. 为了判断端点的敛散性, 输入 Simplifyan/.x-(49/16)则输出右端点处幂级数的一般项为因此,在端点处,级数发散. 再输入 Simplifyan/.x-(47/16)则输出左端点处幂级数的一般项为因此,在端点处, 级数收敛. 也可以在收敛域内求得这个级数的和函数. 输入 Sum4(2n)*(x-3)n/(n+1),n,0,Infinity则输出 函数的幂级数展开 例1.6 (教材 例1.5) 求的6阶麦克劳林展开式. 输入 SeriesCosx,x,0,6则输出 注:这是带皮亚诺余项的麦克劳林展开式. 例1.6 (教材 例1.6) 求在处的6阶泰勒展开式.输入 SeriesLogx,x,1,6则输出例1.7 (教材 例1.7) 求的5阶泰勒展开式.输入serl=SeriesArcTanx,x,0,5;Poly=Normalserl则输出的近似多项式 通过作图把和它的近似多项式进行比较. 输入PlotEvaluateArcTanx,Poly,x,-3/2,3/2,PlotStyle-Dashing0.01,GrayLevel0,AspectRatio-l则输出所作图形（图1.7）, 图中虚线为函数,实线为它的近似多项式.图1.7 例1.9 求在处的8阶泰勒展开, 并通过作图比较函数和它的近似多项式. 输入Clearf;fx_=Exp-(x-1)2*(x+1)2;poly2=NormalSeriesfx,x,1,8PlotEvaluatefx,poly2,x,-1.75,1.75,PlotRange-2,3/2,PlotStyle-Dashing0.01,GrayLevel0则得到近似多项式和它们的图1.8. 图1.8例1.10 求函数在处的阶泰勒展开, 通过作图比较函数和它的近似多项式, 并形成动画进一步观察. 因为 所以输入DoPlotSum(-1)j*x(2j+1)/(2j+1)!,j,0,k,Sinx,x,-40,40,PlotStyle-RGBColor1,0,0,RGBColor0,0,1,k,1,45则输出为的3阶和91阶泰勒展开的图形. 选中其中一幅图形,双击后形成动画. 图1.9是最后一幅图.图1.9例1.11 利用幂级数展开式计算(精确到).因为根据在处的展开式有故前项部分和为输入命令sn_=3(1-1/(5*34)-SumProduct5i-1,i,1,k-1/(5k k!3(4k),k,2,n-1);rn_=Product5i-1,i,1,n-1/5n/n!3(4n-5)/80;delta=10(-10);n0=100;DoPrintn=,n,sn=,Nsn,20;Ifrndelta,Break;Ifn=n0,Printfailed,n,n0则输出结果为 傅里叶级数例1.12 (教材 例1.8) 设是以为周期的周期函数,它在的表达式是 将展开成傅里叶级数. 输入Clearg;gx_:=-1/;-Pi=x0gx_:=1/;0=xPigx_:=gx-2Pi/;PiRGBColor0,1,0; 则输出的图形 (图1.10).图1.10因为是奇函数, 所以它的傅里叶展开式中只含正弦项. 输入b2n_:=b2n=2 Integrate1*Sinn*x,x,0,Pi/Pi;fourier2n_,x_:=Sumb2k*Sink*x,k,1,n;tun_:=Plotgx,Evaluatefourier2n,x, x,-Pi,5 Pi,PlotStyle-RGBColor0,1,0,RGBColor1,0.3,0.5,DisplayFunction-Identity; (*tun是以n为参数的作图命令*)tu2=Tabletun,n,1,30,5; (*tu2是用Table命令作出的6个图形的集合*)toshow=Partitiontu2,2; (*Partition是对集合tu2作分割, 2为分割的参数*)ShowGraphicsArraytoshow (*GraphicsArray是把图形排列的命令*)则输出6个排列着的图形（图1.11）,每两个图形排成一行. 可以看到越大, 的傅里叶级数的前项和与越接近. 图1.11实验2 微分方程（基础实验） 实验目的 理解常微分方程解的概念以及积分曲线和方向场的概念,掌握利用Mathematica求微分方程及方程组解的常用命令和方法. 用Dsolve命令求解微分方程例2.1 (教材 例2.1) 求微分方程 的通解. 输入Clearx,y;DSolvey x+2x*yx=x*Exp-x2,yx,x或 DSolveDyx,x+2x*yx=x*Exp-x2,yx,x则输出微分方程的通解: 其中C1是任意常数. 例2.2 (教材 例2.2) 求微分方程在初始条件下的特解. 输入Clearx,y;DSolvex*y x+yx-Expx=0,y1=2 E,yx,x则输出所求特解: 例2.3 (教材 例2.3) 求微分方程的通解. 输入DSolvey x-2y x+5yx=Expx*Cos2 x,yx,x/Simplify则输出所求通解: 例2.4 (教材 例2.4) 求解微分方程, 并作出其积分曲线.输入g1=TablePlotEx+x3/3+c1+x*c2,x,-5,5,DisplayFunction-Identity,c1,-10,10,5,c2,-5,5,5;Showg1,DisplayFunction-$DisplayFunction;则输出积分曲线的图形（图2.3）.图2.3例2.5 (教材 例2.5) 求微分方程组在初始条件下的特解. 输入Clearx,y,t;DSolvex t+xt+2 yt=Expt, yt -xt- yt=0,x0=1,y0=0,xt,yt,t则输出所求特解: 例2.6 验证是微分方程的通解.输入命令GraphicsPlotFieldn,n,-3,3,x,-3,3;g2=PlotVectorField1,x2/(y4-2),x,-3,3,y,-3,3,Frame-True,ScaleFunction-(1&),ScaleFactor-0.16,HeadLength-0.01,PlotPoints-20,25;g=Showg2,g1,Axes-None,Frame-True;ShowGraphicsArrayg1,g2,g;则分别输出积分曲线如图2.4(a), 微分方程的方向场如图2.4(b). 以及在同一坐标系中画出积分曲线和方向场的图形如下图2.4 (c).(a) (b) (c)图2.4从图2.4(c)中可以看出微分方程的积分曲线与方向场的箭头方向吻合, 且当时, 无论初始条件是什么, 所有的解都趋向于一条直线方程.例2.7 (教材 例2.6) 求解微分方程并作出积分曲线.输入-1,1,-2,2,PlotStyle-RGBColor1,0,0,DisplayFunction-Identity;g2=PlotVectorField1,-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;Showg1,g2,Axes-None,Frame-True,DisplayFunction-$DisplayFunction;则输出积分曲线的图形（图2.5）.图2.5例2.8 求解微分方程并作出其积分曲线.输入命令-3,3,-3,3,PlotStyle-RGBColor1,0,0,DisplayFunction-Identity;gg2=PlotEvaluatet2,x,-3,3,PlotRange-3,3,-3,3,PlotStyle-RGBColor1,0,0,DisplayFunction-Identity;g1=ContourPloty-x3/3-x*(-2+y2),x,-3,3,y,-3,3,PlotRange-3,3,Contours-7,ContourShading-False,PlotPoints-50,DisplayFunction-Identity;g2=PlotVectorField1,(x2+y2-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;Showg1,g2,Axes-None,Frame-True,DisplayFunction-$DisplayFunction;Showgg1,gg2,g2,Axes-None,Frame-True,DisplayFunction-$DisplayFunction;则输出微分方程的向量场与积分曲线, 并输出等值线的图2.6. 图2.6用NDSolve命令求微积分方程的近似解 例2.9 (教材 例2.7) 求初值问题:在区间1.2,4上的近似解并作图. 输入fl=NDSolve(1+x*yx)*yx+(1-x*yx)*yx=0,y1.2=1,y,x,1.2,4则输出为数值近似解(插值函数)的形式: y-InterpolatingFunction1.2,4.,用Plot命令可以把它的图形画出来.不过还需要先使用强制求值命令Evalu-ate, 输入 PlotEvaluateyx/.fl,x,1.2,4则输出近似解的图形（图2.7）.图2.7如果要求区间1.2,4内某一点的函数的近似值, 例如,只要输入 y1.8/.fl则输出所求结果3.8341例2.10 (教材 例2.8) 求范德波尔(Van der Pel)方程在区间0,20上的近似解. 输入Clearx,y;NDSolveyx+(yx2-1)*yx+yx=0,y0=0,y0=-0.5,y,x,0,20;PlotEvaluateyx/.%,x,0,20可以观察到近似解的图形（图2.8）.图2.8的数值解, 并作出数值解的图形.输入命令All,DisplayFunction-Identity;g2=PlotVectorField1,(x2*y*Sinx-1)/x,x,1,4,y,-2,9,Frame-True,ScaleFunction-(1&),ScaleFactor-0.16,HeadLength-0.01,PlotPoints-20,25,DisplayFunction-Identity;g=Showg1,g2,Axes-None,Frame-True;ShowGraphicsArrayg1,g,DisplayFunction-$DisplayFunction;则输出所给微分方程的数值解及数值解的图2.9.图2.9例2.11 (教材 例2.9) 求出初值问题的数值解, 并作出数值解的图形.输入NDSolveyx+Sinx2*yx+yx=Cosx2,y0=1,y0=0,yx,x,0,10PlotEvaluateyx/.%,x,0,10;则输出所求微分方程的数值解及数值解的图形（图2.10）.图2.10例2.12 (教材 例2.10) 洛伦兹(Lorenz)方程组是由三个一阶微分方程组成的方程组.这三个方程看似简单, 也没有包含复杂的函数, 但它的解却很有趣和耐人寻味. 试求解洛伦兹方程组并画出解曲线的图形.输入Cleareq,x,y,zeq=Sequencext=16*yt-16*xt,yt=-xt*zt-yt+45xt,zt=xt*yt-4zt;sol1=NDSolveeq,x0=12,y0=4,z0=0,xt,yt,zt,t,0,16,MaxSteps-10000;g1=ParametricPlot3DEvaluatext,yt,zt/.sol1,t,0,16,PlotPoints-14400,Boxed-False,Axes-None;则输出所求数值解的图形（图2.11(a)）. 从图中可以看出洛伦兹微分方程组具有一个奇异吸引子, 这个吸引子紧紧地把解的图形“吸”在一起. 有趣的是, 无论把解的曲线画得多长, 这些曲线也不相交. (a) (b)图2.11改变初值为输入sol2=NDSolveeq,x0=6,y0=-10,z0=10,xt,yt,zt,t,0,24,MaxSteps-10000;g2=ParametricPlot3DEvaluatext,yt,zt/.sol2,t,0,24,PlotPoints-14400,Boxed-False,Axes-None;ShowGraphicsArrayg1,g2;则输出所求数值解的图形（图2.11(b)）. 从图中可以看出奇异吸引子又出现了, 它把解“吸”在某个区域内, 使得所有的解好象是有规则地依某种模式缠绕.实验3 抛射体的运动（续）（综合实验）实验目的 通过微分方程建模和Mathematica软件,在项目一实验5的基础上,进一步研究在考虑空气阻力的情况下抛射体的运动.问题 根据侦察,发现离我军大炮阵地水平距离10km的前方有一敌军的坦克群正以每小时50km向我军阵地驶来,现欲发射炮弹摧毁敌军坦克群. 为在最短时间内有效摧毁敌军坦克,要求每门大炮都能进行精射击,这样问题就可简化为单门大炮对移动坦克的精确射击问题. 假设炮弹发射速度可控制在0.2km/s至0.5km/s之间,问应选择怎样的炮弹发射速度和怎样的发射角度可以最有效摧毁敌军坦克. 说明 本节我们研究受到重力和空气阻力约束的抛射体的射程. 用记抛射体的位置, 其中x轴是运动的水平方向, y轴是垂直方向. 通过在的约束下最大化x, 可以计算出使抛射体的射程最大的发射角. 假设时抛射体(炮弹)在原点(0,0)以与水平线夹角为初始速度为发射出去. 它受到的空气阻力为 (3.1)重力为 (3.2)在推导和所满足的微分方程之前, 补充一点说明:虽然我们将位置变量仅写作t的函数,但实际上位置变量还依赖于几个其它的变量或参数. 特别是,x和y也依赖于发射角、阻力系数k、质量m及重力加速度g等.为了推导x和y的方程, 按照牛顿定律并结合重力的公式(3.2)和空气阻力的公式(3.1), 得到微分方程 (3.3) (3.4)根据前面所述假设知, 满足下列初始条件, (3.5)先求解,由方程(3.3),令可将其化为一阶微分方程易求出其通解 由 得,所以从通过积分得到x, 即由 得 所以 (3.6)类似地,可从方程(3.4)解出y. 令 方程化为一阶微分方程, 两端除以m,得再在上述方程两端乘以积分因子得即 两端积分得 所以 利用初始条件确定其中的常数C后, 积分v得到y,再次利用初始条件确定任意常数后,则得到 (3.7)下面我们利用公式(3.6)与(3.7)来描绘炮弹运行的典型图形.假定炮弹发射的初速度为0.25km/s, 发射角为, 输入Cleara,t,x,y,g,m,kxv_,a_,t_:=(m/k)*v*Cosa Pi/180*(1-Exp-(k/m)*t)yv_,a_,t_:=(g*m/k)(m/k)-t-(m/k)*Exp-(k/m)*t)+(m/k)*v*Sina Pi/180*(1-Exp-(k/m)*t)g=9.8;m=5.0;k=0.01;炮弹飞行的时间由炮弹落地时的条件所确定. 输入 FindRooty350,55,t=0,t,50则输出炮弹飞行的时间 t-57.4124当发射角时, 输入x350,55, 57.4124/N则输出炮弹的最大射程为 10888.5 现在我们可以画出炮弹运行的典型轨迹了. 输入 ParametricPlotx350,55,t,y350,55,t,t,0,57.4124,PlotRange-0,11000,AxesLabel-x,y则输出图3.1. 图3.1实验报告在上述假设下,进一步研究下列问题: (1) 选择一个初始速度和发射角,利用Mathematica画出炮弹运行的典型轨迹.(2) 假定坦克在大炮前方10km处静止不动,炮弹发射的初速度为0.32km/s,应选择什么样的发射角才能击中坦克?画出炮弹运行的几个轨迹图,通过实验数据和图形来说明你的结论的合理性.(3) 假定坦克在大炮前方10km处静止不动,探索降低或调高炮弹发射的初速度的情况下,应如何选择炮弹的发射角?从上述讨论中总结出最合理有效的发射速度和发射角. (4) 在上题结论的基础上,继续探索,假定坦克在大炮前方10km处以每小时50km向大炮方向前进,此时应如何制定迅速摧毁敌军坦克的方案? 注:在研究过程中,还要包括适当改变阻力系数k与炮弹的质量m所带来的变化.实验4 蹦极跳运动（综合实验）实验目的 利用Mathematica软件,通过微分方程建模,研究蹦极跳运动.问题 在不考虑空气阻力和考虑空气阻力等多种情况下,研究蹦极跳运动中,蹦极者与蹦极绳设计之间的各种关系.说明 蹦极绳相当于一根粗橡皮筋或有弹性的绳子. 当受到张力使之超过其自然长度,绳子会产生一个线性回复力, 即绳子会产生一个力使它恢复到自然长度, 而这个力的大小与它被拉伸的长度成正比. 在一次完美的蹦极跳过程中, 蹦极者爬上一座高桥或高的建筑物, 把绳的一头系在自己身上, 另一头系在一个固定物体如桥栏杆上, 当他跳离桥时, 激动人心的时刻就到来了. 这里要分析的是蹦极者从跳出那一瞬间起他的运动规律.首先要建立坐标系. 假设蹦极者的运动轨迹是垂直的, 因此我们只要用一个坐标来确定他在时刻t的位置. 设y是垂直坐标轴, 单位为英尺, 正向朝下, 选择为桥平面, 时间t的单位为秒, 蹦极者跳出的瞬间为 则表示t时刻蹦极者的位置. 下面我们要求出的表达式.由牛顿第二定律, 物体的质量乘以加速度等于物体所受的力. 我们假设蹦极者所受的力只有重力、空气阻力和蹦极绳产生的回复力. 当然, 直到蹦极者降落的距离大于蹦极绳的自然长度时, 蹦极绳才会产生回复力. 为简单起见, 假设空气阻力的大小与速度成正比, 比例系数为1, 蹦极绳回复力的比例系数为0.4. 这些假设是合理的, 所得到的数学结果与研究所做的蹦极实验非常吻合. 重力加速度现在我们来考虑一次具体的蹦极跳. 假设绳的自然长度为 蹦极者的体重为160lb,则他的质量为斯. 在他到达绳的自然长度(即前, 蹦极者的坠落满足下列初值问题: 利用Mathematica求解上述问题. 输入g=32; m=5; L=200; v1t_,y1t_=vt,yt/.DSolvevt=-g-vt/m,yt=vt,v0=0,y0=0,v,y,t则输出 蹦极者坠落L英尺所用的时间为t1=t/.FindRooty1t=-L,t,24.00609现在我们需要找到当蹦极绳产生回复力后的运动初始条件. 当时, 蹦极者的坠落满足方程初始条件为解初值问题:v2t_,y2t_=vt,yt/.DSolvevt=-g-vt/m-0.4*(L+yt)/m,yt=vt,vt1=v1t1,yt1=-L,v,y,t则输出下列结果 这个解是用复指数函数来表示的.现在蹦极者的位置由命令bungeeyt_=IftAll则输出位置-时间图形(图4.1)图4.1从上图可以看出, 蹦极者在大约13s内由桥面坠落770ft, 然后弹回到桥面下550ft, 上下振动几次, 最终降落到桥面下大约600ft处.实验报告1.在上述问题中(求出需要多长时间蹦极者才能到达他运动轨迹上的最低点, 他能下降到桥面下多少英尺?2.用图描述一个体重为195lb, 用200ft长绳子的蹦极者的坠落. 在绳子对他产生力之前,他能做多长时间的“自由”降落?3.假设你有一根300ft长的蹦极索, 在一组坐标轴上画出你所在实验组的全体成员的运动轨迹草图.4.一个55岁, 体重185lb的蹦极者, 用一根250ft长的蹦极索. 在降落过程中, 他达到的最大速度是多少? 当他最终停止运动时, 他被挂在桥面下多少英尺?5.用不同的空气阻力系数和蹦极索常数做实验, 确定一组合理的参数, 使得在这组参数下, 一个160lb的蹦极者可以回弹到蹦极索的自然长度以上.6.科罗拉多的皇家乔治桥(它跨越皇家乔治峡谷)距谷底1053ft, 一个175lb的蹦极者希望能正好碰到谷底, 则他应使用多长的绳子?7.假如上题中的蹦极者体重增加10lb, 再用同样长的绳子从皇家乔治桥上跳下, 则当他撞到乔治峡谷谷底时, 他的坠落速度是多少?