Mathematica——常微分方程拉氏变换与级数试验

1、国3.5常微分方程、拉氏变换与级数实验学习目标1 .会用Mathematica求解微分方程（组）；2 .能用Mathematica求微分方程（组）的数值解；3 .会利用Mathematica进行拉氏变换与逆变换；4 .能进行幕级数和傅里叶级数的展开。一、常微分方程（组）Mathematica能求常微分方程（组）的准确解，能求解的类型大致覆盖了人工求解的范围,功能很强。但不如人灵活（例如在隐函数和隐方程的处理方面），输出的结果与教材上的答案可能在形式上不同。另外，Mathematica求数值解也很方便，且有利于作出解的图形。在本节中，使用Laplace变换解常微分方程（组）的例子也是十分成功的，

2、过去敬而远之的方法如今可以轻而易举的实现了。求准确解的函数调用格式如下:求方程eqn的通解y（x）,DSolveeqn,yx,x其中自变量是x。DSolveeqn,yx0=y0,yx,x求满足初始条件y（xo）=yo的特解y（x）。DSolveeqn1,eqn2,yix,y2x,x求方程组的通解。DSolveequ1,，yixo=y10,yix,y2x,x求方程组的特解。说明：应当特别注意，方程及各项参数的表述方式很严格，容易出现输入错误。微分方程的表示法只有通过例题才能说清楚。例1解下列常微分方程（组）：5(1)八言+2,(2)y'=1y2(xx3)y'(3)fFy=z,Z=

3、_yrry=zz（0）=1的特解。(4) 3y的通解及潴足初始条件y(0)=0,4=-y解：In1:=DSolvey'x=2yx/(x+1)+(x+1)A(5/2),yx,x0叩My凶号7/2叫In2:=DSolvey'x=(1+yxA2)/(x+xA3)yx),yx,x1-1-2c1x1112x1-1一2c1Out2=yx>rx,yx>112xIn3:=DSolvey'x=zx,z'x=-yx,yx,zx,xOut3=yx一C1Cosx+C2Sinx,zx-C2Cosx-C1SinxIn4:=DSolvey'x=zx,z'x=-yx

4、,y0=0,z0=1,yx,zx,xOut4=yx一Sinx,zx-Cosx提示：认真观察上例，可以从中学习输入格式，未知函数总带有自变量，等号用连续键入两个等号表示，这两点由于不习惯会出错！导数符号用键盘上的撇号，连续两撇表示二阶导数，这与习惯相同。自变量、未知量、初始值的表示法与普通变量相同。说明：输出结果总是尽量用显式解表出，有时反而会使表达式变得复杂，这与教科书的习惯不同。当求显式解遇到问题时，会给出提示。通解中的任意常数用C1,C2,表示。例2求解下列微分方程：(1) y3y3yy=(x-5)e解：In1:=DSolveyx+3yx+3y'(2) x2+(y)2=1,(3)J

5、7=xy。x+yx=(x-5)Exp-x,yx,x222、Out1=yxTe/x2-5x+e'x34e*C1e&xC2ex2C3,34In2：=Simplify%Out2=yx>e(-20x3x424C124xC224x2C3)24In3:=DSolvexA2+y，xA2=1,yx,xOut3=y凶t=-+3，12ArcSinxyx>-x.1-x-2C1In4:=DSolveSqrty'x=xyx,yx,x一一3一Out4=yxx-C1说明：由以上可以看出对方程的类型并无限制，但是输出的答案未必符合习惯，例如第一个方程的答案需要化简，有时即使化简后也未必与教

6、材上的答案一致。例3求微分方程包十2xy=xe'的通解。dx解：In1:=DSolveyx+2xyx=xEa(-xA2),yx,x122。四网yxLfx+及川这就是所给微分方程的通解。式中的C1是通解中的任意常数。上述命令也可以输入为：DSolveDyx+2xyx=xEA（-xA2）,yx,x。例4求微分方程xy+y-ex=0在初始条件y|x=i=2e下的特解。一'一解：In1:=DSolvex*yx+yx-EAx=0,y1=2E,yx,xe-exOut1=yx-二、常微分方程（组）的数值解函数NDSolve用于求给定初值条件或边界条件的常微分方程（组）的近似解，其调用格式如下

7、：NDSolveeqns,y1，y2,x,xmin,xmax求常微分方程（组）的近似解。其中微分方程和初值条件的表示法如同DSolve,未知函数仍有带自变量和不带自变量两种形式，通常使用后一种更方便。初值点xo可以取在区间xmin,xmax上的任何一点处，得到插值函数InterpolatingFunctiondomain,table类型的近似解，近似解的定义域domain一般为domain,table,也有可能缩小。例5求常微分方程y'=x2+y2,满足初始条件y（0）=0的数值解。解：In1:=s1=NDSolvey'x=xA2+yxA2,y0=0,y,x,-2,2Out1=

8、y一InterpolatingFunction-2.,2.,<>In2:=y=y/.s11Out2=InterpolatingFunction-2.,2.,<>In3:=Plotyx,x,-2,2,AspectRatio-Automatic,PlotRange-1.5,1.5图13-43微分方程的解曲线Out3=-Graphics-上例中包含许多值得学习的实用内容，其中第二项参数使用y而不是yx,比用yx好。如果求解区间改为x,-3,3,就会出现警告提示，实际得不到-3,3上的解。Out1表明返回的解放在一个表中，不便使用，实际的解就是插信函数:Interpolatin

9、gFunction-2.,2.,<>In2的结果是用y表示解函数的名字，因此In3顺利画出解曲线如图13-43所示例6求常微分方程组：13x=yxx3y-x满足初始条件x(0)=0,y(0)=1的数值解。解：In1:=s1=NDSolvex't=yt-(xtA3/3-xt),y't=-xt,x0=0,y0=1,x,y,t,-15,15Out1=x一InterpolatingFunction-15.,15.,<>,yfInterpolatingFunction-15.,15.,<>In2：=x=x/.s11,1y=y/.s11,2Out2=In

10、terpolatingFunction-15.,15.,<>Out3=InterpolatingFunction-15.,15.,<>In4：=ParametricPlotxt,yt,t,-15,15,AspectRatio-AutomaticOut3=-Graphics-说明：上例是求一个著名方程组的近似解，其中In2也可以改用一个赋值式x,y=x,y/.Flattens1,一次得到两个函数。通过求数值解容易得到它的相图，In4绘制了解的相轨线如图13-44所示，图中表明原点是奇点，极限环的形状也已经得到。为了应付复杂的情况，需要设置可选参数:WorkingPreci

11、sion参见数值积分部分的介绍。AccuracyGoal计算结果的绝对误差。PrecisionGoal计算结果的相对误差。最大步数。MaxStepsStartingStepSize初始步长。以上可选参数的默认值都为Automatic,其中AccuracyGoal和PrecisionGoal的默认值比WorkingPrecision小10,当解趋于0时应将AccuracyGoal取成Infinity。对于常微分方程，最大步长默认值为1000。这个函数也可以解偏微分方程，最大步长默认值为200。例7解下列微分方程(组)：I .(1) y'=i,潴足初始条件y(0)=1的特解；4yx=3x3

12、y(2)y'=xz+26.5x-y,满足初始条件x(0)=z(0)=0,y(0)=1的特解。z=xy-z解：In1：=NDSolvey'x=I/4yx,y0=1,y,x,1,AccuracyGoal-20,PrecisionGoaH20,WorkingPrecision-25Out1=yTnterpolatingFunction0,1.000000000000000000000000000,<>In2:=y1/.%Out2=0.968912424710644784118519+0.2474039592545229296234109iIn3:=NDSolvex'

13、;t=-3(xt-yt),y't=-xtzt+36.5xt-yt,z't=xtyt-zt,x0=z0=0,y0=1,x,y,z,t,0,20,MaxStepa3000Out3=x一InterpolatingFunction0.,20.,<>,yfInterpolatingFunction0.,20.,<>,zfInterpolatingFunction0.,20.,<>,In4：=ParametricPlot3DEvaluatext,yt,zt/.%,t,0,20,PlotPoints-100010图13-453维相轨线Out3=-Graph

14、ics3D-说明：以上范例中In1取高精度，而且是复系数方程。In2是求解在x=1时的近似值，II 1i求精确解能得到准确值e4,读者可以求e4的近似值与Out2的结果比较，验证近似解的精确度确实很高。In3在求解时增大步数，成功地得到了由In4绘制的如图13-45所示的解的相轨线。In4所示的绘图语句与前面例子中的不同，现在只要会模仿使用它们就行了，要想弄清原理请参阅相关Mathematica书籍。三、拉氏变换Mathematica可以进行拉普拉斯变换，其变换使用的函数调用格式如下：LaplaceTransformf,t,s求函数f(t)的Laplace变换，返回自变量为s的函数。Inver

15、seLaplaceTransformF,s,t求函数F(s)的Laplace逆变换，返回自变量为t的函数。其中函数f(t)和F(s)也可以是函数表，这样可一次变换多个函数。例8求函数t4和etsint的拉氏变换。解：In1：=LaplaceTransformtA4,t,s-24Out1=-5sIn2：=LaplaceTransformExptSint,t,s1Out2=-1r2-2ssIn3:=InverseLaplaceTransform%1,s,t_4Out3=t4In4：=InverseLaplaceTransform%2,s,tOut4=-1iie(14)t(-1e2it)2In5:=

16、FullSimplify%Out5=etSint例9求函数f(t)=t3eat的拉氏变换。解：In1：=LaplaceTransformtA3Expat,t,sOut1=64(-as)以上只是直接进行拉氏变换和逆变换的例子。以下使用拉氏变换解常微分方程，解法原理见本书理论篇，这里完全实现了计算机求解。例10用拉氏变换解微分方程:x“+3x+3x'+x=1满足条件x(0)=x'(0)=x(0)=0的解。解：In1：=f1=LaplaceTransformx'"t+3x"t+3x't+xt,t,sOut1=LaplaceTransformxt,t

17、,s+s3LaplaceTransformxt,t,s+3(sLaplaceTransformxt,t,s-x0)-s2x0+3(s2LaplaceTransformxt,t,s-sx0-x'0)-sx'0-x0In2：=s1=LaplaceTransform1,t,s-1Out2=sIn3:=x0=x，0=x0=0;Solvef1=s1,LaplaceTransformxt,t,s1Out4=LaplaceTransformxt,t,s3s(1s)s,t1In5:=InverseLaplaceTransform3,s(1s)1j_2Out5=1-ex(22tt2)说明：上例中

18、的LaplaceTransformxt,t,s就是教材中的X(s),In3解出X(s),其余过程与教科书完全相同。现在可以将一切计算留给计算机，学生只要弄清解法原理及过程。技巧：充分利用复制、粘贴功能，可以加快输入速度，避免键入错误。上例中In5就可以从Out4中将表达式复制过来。例11求微分方程组：父-2x'-y'+2y=0、x'+y.-2x=-2e上满足条件x(0)=3,x(0)=2,y(0)=0的特解。解：In1：=f1=LaplaceTransformxt-2x't-y't+2yt,x't+y't-2xt,t,s;In2：=s1=

19、LaplaceTransform0,-2Exp-t,t,s;In3:=x0=3;x'0=2；y0=0;Solvef1=s1,LaplaceTransformxt,t,s,LaplaceTransformyt,t,s;In5:=InverseLaplaceTransformFlattenLaplaceTransformxt,t,s,LaplaceTransformyt,t,s/.%,s,tOut5=5-e-t-3et+2e2：e-t(-1+et)2(1+2et)In6:=Simplify%Out6=5-e-t-3et+2e2t,e-t-3et+2e2t说明：在上例中，不显示任何中间结果，

20、语句比较简练。其中，In1和In2分别对方程组的左边和右边进行拉氏变换，In3解出X(s)和Y(s)。In5比较难懂，可以参看前面的例题，这里是从Out3中自动将解X(s)和Y(s)提取出来，再进行拉氏逆变换。Out5是x(t),y(t),Out6将答案化简。本例已经将求解过程一般化，只需改变方程组和初值的数据，就可以解其它方程组了。四、级数1 .求和与求积求有限或无穷和、积的函数是：imaxSumf,i,imin,imax求工f(i),其中imin可以是-00,i4minimax可以是00(即+00),但是必须满足imin&imax。基本输入模板中也有求和专用的符号，使用模板输入更方

21、便。求多重和，也可以使用基本输入模imax求口f(i),基本输入模板中也有i=imin求多重积，也可以使用基本输入Sumf,i,imin,imax,j,jmin,jmax,板连续多次输入求和符号得到。Productf,i,imin,imax求积符号Productf,i,imin,imax,j,jmin,jmax,模板连续多次输入求积符号得到。二1,、二1£，(4)口ekk=1kk=1例12求下列级数的和与积：nc二1(1)£k2,(2)£(3)k1kWk解：In1：=SumkA2,k,1,n1Out1=n(1n)(12n)6coIn2:=Z1/kA2k1Out2=

22、odIn3:=Z1/kk1Sum：:div:Sumdoesnotconverge.cdIn4:=nExp1/kA2k1Out4=e手说明：上例中第三个级数发散，Mathematica给出提示，并在不能给出结果时将输入的式子作为输出。NSum和NProduct得到数值解。2 .将函数展开为幕级数将函数展开为幕级数的函数调用格式如下:Seriesf,x,xo,n为止。将函数f(x)在X0处展成幕级数直到n次项Seriesfx,xo,n,y,y0,m例13展开下列函数为幕级数：将函数f(x,y)先对y后对x展开。sinx(1)y=tgx,(2)y=(3) y=f(x),/、xy(4) y=e。解：I

23、n1：=SeriesTanx,x,0,9-x32x517x7Out1=x31531562x92835ox10In2：=SeriesSinx/x,x,0,9Out2=1-+_+6120504036288010oxIn3:=Seriesfx,x,1,71.O1(Q)QOut3=f1f1(x-1)-f1(x-1)-f()1(x-1)214f/)120f(5)1(xF5+j)1(x-1)6f(7)1(x-1)7ox.187205040In4：=SeriesExpxy,x,0,3,y,0,22Out4=1(yoy3)x-yoy3x2oy3x3ox42说明：上例中In3表明也可以展开抽象的函数。对已经展开

24、的幕级数进行操作的两个函数是：Normalexpr将幕级数expr去掉余项转换成多项式SeriesCoefficientexpr,n找出幕级数expr的n次项系数。例14将y=arcsinx展开为幕级数，只取前9项并去掉余项。解：In1：=SeriesArcSinx,x,0,93579Out1=10oxx3x5x35xx-6401121152In2：=Normal%Out2=3o5x3xx6405x7+11235x91152In3:=SeriesCoefficient%1,5一3Out3=403 .傅里叶级数求傅里叶级数就是求出傅里叶系数，傅里叶系数是一个积分表达式，所以利用积分函数Integrate就可以实现。例如，设周期矩形脉冲信号的脉冲宽度为r,脉冲幅度为E,周期为T,这种信号在一1t悍a|t|