Mathematica——常微分方程拉氏變換與級數(shù)試驗

1、國3.5常微分方程、拉氏變換與級數(shù)實驗學習目標1 .會用Mathematica求解微分方程（組）；2 .能用Mathematica求微分方程（組）的數(shù)值解；3 .會利用Mathematica進行拉氏變換與逆變換；4 .能進行幕級數(shù)和傅里葉級數(shù)的展開。一、常微分方程（組）Mathematica能求常微分方程（組）的準確解，能求解的類型大致覆蓋了人工求解的范圍,功能很強。但不如人靈活（例如在隱函數(shù)和隱方程的處理方面），輸出的結(jié)果與教材上的答案可能在形式上不同。另外，Mathematica求數(shù)值解也很方便，且有利于作出解的圖形。在本節(jié)中，使用Laplace變換解常微分方程（組）的例子也是十分成功的，

2、過去敬而遠之的方法如今可以輕而易舉的實現(xiàn)了。求準確解的函數(shù)調(diào)用格式如下:求方程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求方程組的特解。說明：應當特別注意，方程及各項參數(shù)的表述方式很嚴格，容易出現(xiàn)輸入錯誤。微分方程的表示法只有通過例題才能說清楚。例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提示：認真觀察上例，可以從中學習輸入格式，未知函數(shù)總帶有自變量，等號用連續(xù)鍵入兩個等號表示，這兩點由于不習慣會出錯！導數(shù)符號用鍵盤上的撇號，連續(xù)兩撇表示二階導數(shù)，這與習慣相同。自變量、未知量、初始值的表示法與普通變量相同。說明：輸出結(jié)果總是盡量用顯式解表出，有時反而會使表達式變得復雜，這與教科書的習慣不同。當求顯式解遇到問題時，會給出提示。通解中的任意常數(shù)用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。四網(wǎng)yxLfx+及川這就是所給微分方程的通解。式中的C1是通解中的任意常數(shù)。上述命令也可以輸入為：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-二、常微分方程（組）的數(shù)值解函數(shù)NDSolve用于求給定初值條件或邊界條件的常微分方程（組）的近似解，其調(diào)用格式如下

7、：NDSolveeqns,y1，y2,x,xmin,xmax求常微分方程（組）的近似解。其中微分方程和初值條件的表示法如同DSolve,未知函數(shù)仍有帶自變量和不帶自變量兩種形式，通常使用后一種更方便。初值點xo可以取在區(qū)間xmin,xmax上的任何一點處，得到插值函數(shù)InterpolatingFunctiondomain,table類型的近似解，近似解的定義域domain一般為domain,table,也有可能縮小。例5求常微分方程y'=x2+y2,滿足初始條件y（0）=0的數(shù)值解。解：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-上例中包含許多值得學習的實用內(nèi)容，其中第二項參數(shù)使用y而不是yx,比用yx好。如果求解區(qū)間改為x,-3,3,就會出現(xiàn)警告提示，實際得不到-3,3上的解。Out1表明返回的解放在一個表中，不便使用，實際的解就是插信函數(shù):Interpolatin

9、gFunction-2.,2.,<>In2的結(jié)果是用y表示解函數(shù)的名字，因此In3順利畫出解曲線如圖13-43所示例6求常微分方程組：13x=yxx3y-x滿足初始條件x(0)=0,y(0)=1的數(shù)值解。解：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,一次得到兩個函數(shù)。通過求數(shù)值解容易得到它的相圖，In4繪制了解的相軌線如圖13-44所示，圖中表明原點是奇點，極限環(huán)的形狀也已經(jīng)得到。為了應付復雜的情況，需要設置可選參數(shù):WorkingPreci

11、sion參見數(shù)值積分部分的介紹。AccuracyGoal計算結(jié)果的絕對誤差。PrecisionGoal計算結(jié)果的相對誤差。最大步數(shù)。MaxStepsStartingStepSize初始步長。以上可選參數(shù)的默認值都為Automatic,其中AccuracyGoal和PrecisionGoal的默認值比WorkingPrecision小10,當解趨于0時應將AccuracyGoal取成Infinity。對于常微分方程，最大步長默認值為1000。這個函數(shù)也可以解偏微分方程，最大步長默認值為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取高精度，而且是復系數(shù)方程。In2是求解在x=1時的近似值，II 1i求精確解能得到準確值e4,讀者可以求e4的近似值與Out2的結(jié)果比較，驗證近似解的精確度確實很高。In3在求解時增大步數(shù)，成功地得到了由In4繪制的如圖13-45所示的解的相軌線。In4所示的繪圖語句與前面例子中的不同，現(xiàn)在只要會模仿使用它們就行了，要想弄清原理請參閱相關Mathematica書籍。三、拉氏變換Mathematica可以進行拉普拉斯變換，其變換使用的函數(shù)調(diào)用格式如下：LaplaceTransformf,t,s求函數(shù)f(t)的Laplace變換，返回自變量為s的函數(shù)。Inver

15、seLaplaceTransformF,s,t求函數(shù)F(s)的Laplace逆變換，返回自變量為t的函數(shù)。其中函數(shù)f(t)和F(s)也可以是函數(shù)表，這樣可一次變換多個函數(shù)。例8求函數(shù)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求函數(shù)f(t)=t3eat的拉氏變換。解：In1：=LaplaceTransformtA3Expat,t,sOut1=64(-as)以上只是直接進行拉氏變換和逆變換的例子。以下使用拉氏變換解常微分方程，解法原理見本書理論篇，這里完全實現(xiàn)了計算機求解。例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),其余過程與教科書完全相同。現(xiàn)在可以將一切計算留給計算機，學生只要弄清解法原理及過程。技巧：充分利用復制、粘貼功能，可以加快輸入速度，避免鍵入錯誤。上例中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說明：在上例中，不顯示任何中間結(jié)果，

20、語句比較簡練。其中，In1和In2分別對方程組的左邊和右邊進行拉氏變換，In3解出X(s)和Y(s)。In5比較難懂，可以參看前面的例題，這里是從Out3中自動將解X(s)和Y(s)提取出來，再進行拉氏逆變換。Out5是x(t),y(t),Out6將答案化簡。本例已經(jīng)將求解過程一般化，只需改變方程組和初值的數(shù)據(jù)，就可以解其它方程組了。四、級數(shù)1 .求和與求積求有限或無窮和、積的函數(shù)是：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,板連續(xù)多次輸入求和符號得到。Productf,i,imin,imax求積符號Productf,i,imin,imax,j,jmin,jmax,模板連續(xù)多次輸入求積符號得到。二1,、二1£，(4)口ekk=1kk=1例12求下列級數(shù)的和與積：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手說明：上例中第三個級數(shù)發(fā)散，Mathematica給出提示，并在不能給出結(jié)果時將輸入的式子作為輸出。NSum和NProduct得到數(shù)值解。2 .將函數(shù)展開為幕級數(shù)將函數(shù)展開為幕級數(shù)的函數(shù)調(diào)用格式如下:Seriesf,x,xo,n為止。將函數(shù)f(x)在X0處展成幕級數(shù)直到n次項Seriesfx,xo,n,y,y0,m例13展開下列函數(shù)為幕級數(shù)：將函數(shù)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表明也可以展開抽象的函數(shù)。對已經(jīng)展開

24、的幕級數(shù)進行操作的兩個函數(shù)是：Normalexpr將幕級數(shù)expr去掉余項轉(zhuǎn)換成多項式SeriesCoefficientexpr,n找出幕級數(shù)expr的n次項系數(shù)。例14將y=arcsinx展開為幕級數(shù)，只取前9項并去掉余項。解：In1：=SeriesArcSinx,x,0,93579Out1=10oxx3x5x35xx-6401121152In2：=Normal%Out2=3o5x3xx6405x7+11235x91152In3:=SeriesCoefficient%1,5一3Out3=403 .傅里葉級數(shù)求傅里葉級數(shù)就是求出傅里葉系數(shù)，傅里葉系數(shù)是一個積分表達式，所以利用積分函數(shù)Integrate就可以實現(xiàn)。例如，設周期矩形脈沖信號的脈沖寬度為r,脈沖幅度為E,周期為T,這種信號在一1t悍a|t|