Mathematica——常微分方程、拉氏變換與級(jí)數(shù)實(shí)驗(yàn).doc

13.5 常微分方程、拉氏變換與級(jí)數(shù)實(shí)驗(yàn)學(xué)習(xí)目標(biāo)1. 會(huì)用Mathematica求解微分方程（組）；2. 能用Mathematica求微分方程（組）的數(shù)值解；3. 會(huì)利用Mathematica進(jìn)行拉氏變換與逆變換；4. 能進(jìn)行冪級(jí)數(shù)和傅里葉級(jí)數(shù)的展開(kāi)。一、 常微分方程（組） Mathematica能求常微分方程（組）的準(zhǔn)確解，能求解的類型大致覆蓋了人工求解的范圍，功能很強(qiáng)。但不如人靈活（例如在隱函數(shù)和隱方程的處理方面），輸出的結(jié)果與教材上的答案可能在形式上不同。另外，Mathematica求數(shù)值解也很方便，且有利于作出解的圖形。在本節(jié)中，使用Laplace變換解常微分方程（組）的例子也是十分成功的，過(guò)去敬而遠(yuǎn)之的方法如今可以輕而易舉的實(shí)現(xiàn)了。 求準(zhǔn)確解的函數(shù)調(diào)用格式如下： DSolveeqn，yx，x 求方程eqn的通解y（x），其中自變量是x。 DSolveeqn，yx0= =y0，yx，x 求滿足初始條件y（x0）= y0的特解y（x）。 DSolveeqn1，eqn2，y1x，y2x，x 求方程組的通解。 DSolveequ1，y1x0= =y10，y1x，y2x，x 求方程組的特解。 說(shuō)明：應(yīng)當(dāng)特別注意，方程及各項(xiàng)參數(shù)的表述方式很嚴(yán)格，容易出現(xiàn)輸入錯(cuò)誤。微分方程的表示法只有通過(guò)例題才能說(shuō)清楚。例1 解下列常微分方程（組）： （1），（2）， （3） ， （4）的通解及滿足初始條件y（0）=0，z（0）=1的特解。 解：In1：=DSolveyx= =2yx/（x+1）+（x+1）（5/2）， yx，x Out1= In2：=DSolveyx= =（1+yx2）/(x+x3）yx)，yx，x Out2=， In3：=DSolveyx= =zx，zx= = -yx， yx，zx，x Out3=yxC1Cosx+ C2Sinx， zxC2Cosx- C1Sinx In4：=DSolveyx= =zx，zx= = -yx，y0= =0，z0= =1， yx，zx，x Out4=yxSinx，zxCosx 提示：認(rèn)真觀察上例，可以從中學(xué)習(xí)輸入格式，未知函數(shù)總帶有自變量，等號(hào)用連續(xù)鍵入兩個(gè)等號(hào)表示，這兩點(diǎn)由于不習(xí)慣會(huì)出錯(cuò)！導(dǎo)數(shù)符號(hào)用鍵盤(pán)上的撇號(hào)，連續(xù)兩撇表示二階導(dǎo)數(shù)，這與習(xí)慣相同。自變量、未知量、初始值的表示法與普通變量相同。 說(shuō)明：輸出結(jié)果總是盡量用顯式解表出，有時(shí)反而會(huì)使表達(dá)式變得復(fù)雜，這與教科書(shū)的習(xí)慣不同。當(dāng)求顯式解遇到問(wèn)題時(shí)，會(huì)給出提示。通解中的任意常數(shù)用C1，C2，表示。例2 求解下列微分方程： （1），（2），（3）。解：In1：=DSolve+3yx +3yx + yx = =（x - 5）Exp-x， yx，x Out1= In2：=Simplify% Out2= In3：=DSolvex2 + yx2 = = 1，yx，x Out3=， In4：=DSolveSqrtyx = = x yx，yx，x Out4= 說(shuō)明：由以上可以看出對(duì)方程的類型并無(wú)限制，但是輸出的答案未必符合習(xí)慣，例如第一個(gè)方程的答案需要化簡(jiǎn)，有時(shí)即使化簡(jiǎn)后也未必與教材上的答案一致。例3 求微分方程的通解。 解：In1：=DSolveyx+2x yx= = x E（-x2），yx，x Out1=yx 這就是所給微分方程的通解。式中的C1是通解中的任意常數(shù)。 上述命令也可以輸入為：DSolveDyx + 2x yx= =x E（ - x2），yx，x。例4 求微分方程xy+ y - ex = 0在初始條件y|x=1 = 2e下的特解。 解：In1：=DSolvex*yx+yx-Ex= =0，y1= =2E，yx，x Out1= yx二、 常微分方程（組）的數(shù)值解 函數(shù)NDSolve用于求給定初值條件或邊界條件的常微分方程（組）的近似解，其調(diào)用格式如下： NDSolveeqns，y1，y2，x，xmin，xmax 求常微分方程（組）的近似解。 其中微分方程和初值條件的表示法如同DSolve，未知函數(shù)仍有帶自變量和不帶自變量?jī)煞N形式，通常使用后一種更方便。初值點(diǎn)x0可以取在區(qū)間xmin，xmax上的任何一點(diǎn)處，得到插值函數(shù)InterpolatingFunctiondomain，table類型的近似解，近似解的定義域domain一般為domain，table，也有可能縮小。例5 求常微分方程y= x2 + y2，滿足初始條件y（0）= 0的數(shù)值解。解：In1：=s1=NDSolveyx= =x2+yx2，y0= =0， y，x，-2，2 Out1=yInterpolatingFunction-2.，2.， In2：= y=y / . s11 Out2=InterpolatingFunction-2.，2.， In3：=Plotyx，x，-2，2，AspectRatioAutomatic， PlotRange-1.5，1.5圖13-43 微分方程的解曲線 Out3= -Graphics- 上例中包含許多值得學(xué)習(xí)的實(shí)用內(nèi)容，其中第二項(xiàng)參數(shù)使用y而不是yx，比用yx好。如果求解區(qū)間改為x，-3，3，就會(huì)出現(xiàn)警告提示，實(shí)際得不到-3，3上的解。Out1表明返回的解放在一個(gè)表中，不便使用，實(shí)際的解就是插值函數(shù)： InterpolatingFunction-2.，2.， In2的結(jié)果是用y表示解函數(shù)的名字，因此In3順利畫(huà)出解曲線如圖13-43所示。例6 求常微分方程組： 滿足初始條件x（0）=0，y（0）=1的數(shù)值解。 解：In1：=s1=NDSolvext= = yt -（xt3/3 - xt）， yt= = - xt，x0= =0，y0= =1， x，y，t，-15，15 Out1=xInterpolatingFunction-15.，15.， yInterpolatingFunction-15.，15.， In2：= x=x / . s11，1 y=y / . s11，2 Out2=InterpolatingFunction-15.，15.， Out3=InterpolatingFunction-15.，15.， In4：=ParametricPlotxt，yt，t，-15，15， AspectRatioAutomatic圖13-44 解的相軌線 Out3= -Graphics- 說(shuō)明：上例是求一個(gè)著名方程組的近似解，其中In2也可以改用一個(gè)賦值式x，y=x，y / . Flattens1，一次得到兩個(gè)函數(shù)。通過(guò)求數(shù)值解容易得到它的相圖，In4繪制了解的相軌線如圖13-44所示，圖中表明原點(diǎn)是奇點(diǎn)，極限環(huán)的形狀也已經(jīng)得到。 為了應(yīng)付復(fù)雜的情況，需要設(shè)置可選參數(shù)： WorkingPrecision 參見(jiàn)數(shù)值積分部分的介紹。 AccuracyGoal 計(jì)算結(jié)果的絕對(duì)誤差。 PrecisionGoal 計(jì)算結(jié)果的相對(duì)誤差。 MaxSteps 最大步數(shù)。 StartingStepSize 初始步長(zhǎng)。 以上可選參數(shù)的默認(rèn)值都為Automatic，其中AccuracyGoal和PrecisionGoal的默認(rèn)值比WorkingPrecision小10，當(dāng)解趨于0時(shí)應(yīng)將AccuracyGoal取成Infinity。對(duì)于常微分方程，最大步長(zhǎng)默認(rèn)值為1000。這個(gè)函數(shù)也可以解偏微分方程，最大步長(zhǎng)默認(rèn)值為200。例7 解下列微分方程（組）：（1），滿足初始條件y（0）=1的特解；（2），滿足初始條件x（0）=z（0）=0，y（0）=1的特解。 解：In1：=NDSolveyx= =I/4yx，y0= =1，y，x，1， AccuracyGoal20，PrecisionGoal20，WorkingPrecision25 Out1=yInterpolatingFunction 0，1.000000000000000000000000000， In2：=y1 / . % Out2=0.968912424710644784118519 + 0.2474039592545229296234109 In3：=NDSolvext= = -3（xt -yt）， yt = = -xt zt+36.5xt -yt， zt = = xt yt- zt， x0 = = z0 = = 0，y0= =1， x，y，z，t，0，20，MaxSteps3000 Out3=xInterpolatingFunction0.，20.， yInterpolatingFunction0.，20.， zInterpolatingFunction0.，20.， In4：=ParametricPlot3DEvaluatext，yt，zt / . %， t，0，20，PlotPoints 1000圖13-45 3維相軌線 Out3= -Graphics3D-說(shuō)明：以上范例中In1取高精度，而且是復(fù)系數(shù)方程。In2是求解在x=1時(shí)的近似值，求精確解能得到準(zhǔn)確值，讀者可以求的近似值與Out2的結(jié)果比較，驗(yàn)證近似解的精確度確實(shí)很高。In3在求解時(shí)增大步數(shù)，成功地得到了由In4繪制的如圖13-45所示的解的相軌線。In4所示的繪圖語(yǔ)句與前面例子中的不同，現(xiàn)在只要會(huì)模仿使用它們就行了，要想弄清原理請(qǐng)參閱相關(guān)Mathematica書(shū)籍。三、 拉氏變換 Mathematica可以進(jìn)行拉普拉斯變換，其變換使用的函數(shù)調(diào)用格式如下： LaplaceTransformf，t，s 求函數(shù)f（t）的Laplace變換，返回自變量為s的函數(shù)。 InverseLaplaceTransformF，s，t 求函數(shù)F（s）的Laplace逆變換，返回自變量為t的函數(shù)。 其中函數(shù)f（t）和F（s）也可以是函數(shù)表，這樣可一次變換多個(gè)函數(shù)。例8 求函數(shù)t 4和et sint的拉氏變換。 解：In1：=LaplaceTransformt4，t，s Out1= In2：=LaplaceTransformExpt Sint，t，s Out2= In3：=InverseLaplaceTransform%1，s，t Out3=t4 In4：=InverseLaplaceTransform%2，s，t Out4= In5：=FullSimplify% Out5=et Sint例9 求函數(shù)f（t）= t3 eat的拉氏變換。 解：In1：= LaplaceTransformt3 Expa t，t，s Out1= 以上只是直接進(jìn)行拉氏變換和逆變換的例子。以下使用拉氏變換解常微分方程，解法原理見(jiàn)本書(shū)理論篇，這里完全實(shí)現(xiàn)了計(jì)算機(jī)求解。例10 用拉氏變換解微分方程： + 3x+ 3x+ x = 1滿足條件x(0) = x(0) = x(0) = 0的解。解：In1：=f1=LaplaceTransform t +3xt+ 3xt+xt， t，s Out1=LaplaceTransformxt，t，s + s3LaplaceTransformxt，t，s + 3（sLaplaceTransformxt，t，s - x0）- s2x0 + 3（s2LaplaceTransformxt，t，s - s x0 - x0）- s x0- x0 In2：=s1=LaplaceTransform1，t，s Out2= In3：= x0= x0= x0=0； Solvef1= =s1，LaplaceTransformxt，t，s Out4= In5：=InverseLaplaceTransform，s，t Out5= 說(shuō)明：上例中的LaplaceTransformxt，t，s就是教材中的X（s），In3解出X（s），其余過(guò)程與教科書(shū)完全相同?，F(xiàn)在可以將一切計(jì)算留給計(jì)算機(jī)，學(xué)生只要弄清解法原理及過(guò)程。 技巧：充分利用復(fù)制、粘貼功能，可以加快輸入速度，避免鍵入錯(cuò)誤。上例中In5就可以從Out4中將表達(dá)式復(fù)制過(guò)來(lái)。例11 求微分方程組： 滿足條件x（0）=3，x（0）=2，y（0）=0的特解。解：In1：=f1=LaplaceTransform xt - 2xt - yt + 2yt，xt + yt - 2xt， t，s； In2：=s1=LaplaceTransform0，- 2Exp-t，t，s； In3：= x0=3； x0= 2；y0=0； Solvef1= =s1，LaplaceTransformxt，t，s， LaplaceTransformyt，t，s； In5：=InverseLaplaceTransform Flatten LaplaceTransformxt，t，s， LaplaceTransformyt，t，s / . %，s，t Out5=5 - e-t - 3et + 2e2t，e-t（- 1 + et）2（1 + 2 et） In6：=Simplify% Out6= 5 - e-t - 3et + 2e2t，e-t - 3et + 2e2t 說(shuō)明：在上例中，不顯示任何中間結(jié)果，語(yǔ)句比較簡(jiǎn)練。其中，In1和In2分別對(duì)方程組的左邊和右邊進(jìn)行拉氏變換，In3解出X（s）和Y（s）。In5比較難懂，可以參看前面的例題，這里是從Out3中自動(dòng)將解X（s）和Y（s）提取出來(lái)，再進(jìn)行拉氏逆變換。Out5是x（t），y（t），Out6將答案化簡(jiǎn)。本例已經(jīng)將求解過(guò)程一般化，只需改變方程組和初值的數(shù)據(jù)，就可以解其它方程組了。四、 級(jí)數(shù)1 求和與求積 求有限或無(wú)窮和、積的函數(shù)是： Sumf，i，imin，imax 求，其中imin可以是-，imax可以是（即+），但是必須滿足iminimax?；据斎肽０逯幸灿星蠛蛯Ｓ玫姆?hào)，使用模板輸入更方便。 Sumf，i，imin，imax，j，jmin，jmax， 求多重和，也可以使用基本輸入模板連續(xù)多次輸入求和符號(hào)得到。 Productf，i，imin，imax 求，基本輸入模板中也有求積符號(hào)。 Productf，i，imin，imax，j，jmin，jmax， 求多重積 ，也可以使用基本輸入模板連續(xù)多次輸入求積符號(hào)得到。例12 求下列級(jí)數(shù)的和與積： （1），（2） ，（3） ，（4） 。 解：In1：=Sumk2，k，1，n Out1= In2：= Out2= In3：= Sum：div：Sum does not converge. Out3= In4：= Out4= 說(shuō)明：上例中第三個(gè)級(jí)數(shù)發(fā)散，Mathematica給出提示，并在不能給出結(jié)果時(shí)將輸入的式子作為輸出。 NSum和NProduct得到數(shù)值解。2 將函數(shù)展開(kāi)為冪級(jí)數(shù)將函數(shù)展開(kāi)為冪級(jí)數(shù)的函數(shù)調(diào)用格式如下： Seriesf，x，x0，n 將函數(shù)f（x）在x0 處展成冪級(jí)數(shù)直到n次項(xiàng)為止。 Seriesf，x，x0，n，y，y0，m 將函數(shù)f（x，y）先對(duì)y后對(duì)x展開(kāi)。例13 展開(kāi)下列函數(shù)為冪級(jí)數(shù)：（1） y=tgx，（2） ， （3）y = f（x），（4）y = exy。 解：In1：=SeriesTanx，x，0，9 Out1= In2：=SeriesSinx /x，x，0，9 Out2= In3：=Seriesfx，x，1，7 Out3= In4：=SeriesExpx y，x，0，3，y，0，2 Out4= 說(shuō)明：上例中In3表明也可以展開(kāi)抽象的函數(shù)。 對(duì)已經(jīng)展開(kāi)的冪級(jí)數(shù)進(jìn)行操作的兩個(gè)函數(shù)是： Normalexpr 將冪級(jí)數(shù)expr去掉余項(xiàng)轉(zhuǎn)換成多項(xiàng)式。 SeriesCoefficientexpr，n 找出冪級(jí)數(shù)expr的n次項(xiàng)系數(shù)。例14 將y = arcsinx展開(kāi)為冪級(jí)數(shù)，只取前9項(xiàng)并去掉余項(xiàng)。 解：In1：=SeriesArcSinx，x，0，9 Out1= In2：=Normal% Out2= In3：=SeriesCoefficient