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

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

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

3、x x3）y，(3)r fy =z=Z=_yy = z（4）丿'的通解及满足初始条件y （0） =0，z （0） =1的特解。工=-y解：In1: =DSolvey ' x= =2yx/ (x+1) + (x+1) A (5/2)，yx，xf rOut1= "yxT 2(1+x)7/2+(1+x)2c13ln2: =DSolvey ' x= = (1+yxF2 ) /(x+xA3 ) yx)，yx，x In3: =DSolvey ' x= =zx , z' x= = -yx, yx , zx , xOut2= yx >Out3=yx C1

4、Cosx+ C2Sinx,zx C2Cosx- C1Si nxIn4: =DSolvey ' x= =zx , z' x= = -yx , y0= =0 , z0= =1,yx , zx , xOut4=yx Sinx, zx Cosx提示：认真观察上例，可以从中学习输入格式，未知函数总带有自变量，等号用连续键 入两个等号表示，这两点由于不习惯会出错！导数符号用键盘上的撇号，连续两撇表示二阶 导数，这与习惯相同。自变量、未知量、初始值的表示法与普通变量相同。表示。说明：输出结果总是尽量用显式解表出，有时反而会使表达式变得复杂，这与教科书的 习惯不同。当求显式解遇到问题时，会给出

5、提示。通解中的任意常数用C1，C2，例2求解下列微分方程：(1) y 3y 3y y=(x-5)e，(2) x2(y)2 =1，(3) . y 二 xy。解：In1: =DSolve y x +3y x +3y ' x + yx = = (x - 5) Exp-x，yx，xOut1=/ 2、Lx5x23、x_5x + + e x< 2丿123丿yx > ex22-34、-区 xe»C1 exC2 ex2C3 .3 4In2: =Simplify%1Out2= yx e(-20x3 x4 24C1 24xC2 24x2C3) 24In3: =DSolvexA2 +

6、y ' xA2 = = 1，yx，xOut3= yx-Ixd -x ArcSinx C1，2 2泊知1-八。1ln4: =DSolveSqrty ' x = = x yx ，yx，x3Out4= yx 3 -x -C1说明：由以上可以看出对方程的类型并无限制，但是输出的答案未必符合习惯，例如第 一个方程的答案需要化简，有时即使化简后也未必与教材上的答案一致。例3求微分方程 鱼'2xy =xe"的通解。dx解：In1: =DSolvey x+2x yx= = x Ea (帜八2)，yx，x1 2 2 2°ut1=yx-尹 x 宀这就是所给微分方程的通解

7、。式中的C1是通解中的任意常数。上述命令也可以输入为：DSolveDyx + 2x yx= =x EA（ - xA2）, yx , x。例4求微分方程xy + y - ex = 0在初始条件y|x=i = 2e下的特解。解：In1: =DSolvex*y x+yx-EAx= =0 ，y1= =2E，yx，xe +exOut1= yx-x二、常微分方程（组）的数值解函数NDSolve用于求给定初值条件或边界条件的常微分方程（组）的近似解，其调用格 式如下：NDSolveeqns, y 1, y2, x , xmin , xmax求常微分方程（组）的近似解。其中微分方程和初值条件的表示法如同DSo

8、lve，未知函数仍有带自变量和不带自变量两种形式，通常使用后一种更方便。初值点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=y InterpolatingFunction-2. , 2. , <

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

10、tion-2. , 2., < >ln2的结果是用y表示解函数的名字，因此In3顺利画出解曲线如图13-43所示例6 求常微分方程组：13x y x x3y " = x满足初始条件x（0）=0，y（0）=1的数值解。解：In1: =$仁NDSolvex ' t= = yt - （xtA3/3 - xt）, y' t= = - xt , x0= =0 , y0= =1, x , y, t, -15, 15Out1=x f InterpolatingFunction-15. , 15., < >, y InterpolatingFunction-1

11、5. , 15. , < >In2: = x=x / . s11 , 1y=y / . s11 , 2Out2=InterpolatingFunction-15. , 15., < >Out3=InterpolatingFunction-15. , 15., < > ln4: =ParametricPlotxt , yt , t, -15 , 15, AspectRatio Automatic69Out3= -Graphics-说明：上例是求一个著名方程组的近似解，其中 In2也可以改用一个赋值式x , y=x , y / . Flattens1, 一次得到

12、两个函数。通过求数值解容易得到它的相图,In4绘制了解的相轨 线如图13-44所示，图中表明原点是奇点，极限环的形状也已经得到。为了应付复杂的情况，需要设置可选参数:Worki ngPrecisio n参见数值积分部分的介绍AccuracyGoal计算结果的绝对误差Precisi on Goal计算结果的相对误差MaxSteps最大步数。Starti ngStepSize初始步长。以上可选参数的默认值都为 Automatic，其中AccuracyGoal和PrecisionGoal的默认值比 WorkingPrecision小10，当解趋于0时应将AccuracyGoal取成Infinity。

13、对于常微分方程，最 大步长默认值为1000。这个函数也可以解偏微分方程，最大步长默认值为 200。例7解下列微分方程（组）：1（1）y' i，满足初始条件y （0）=1的特解；4yX x = _3x 3y（2）収=xz+26.5x-y，满足初始条件 x（0）=z（0）=0, y（0）=1的特解。z = xy z解：In1: =NDSolvey ' x= =l/4yx ，y0= =1，y，x， 1，AccuracyGoal20， PrecisionGoa220， WorkingPrecision25Out1=y In terpolat ingFun cti on0，1.00000

14、0000000000000000000000, < >In2: =y1 / . %Out2=0.968912424710644784118519 +0.24740395925452292962341091ln3: =NDSolvex ' t= = -3 （xt -yt），y' t = = -xt zt+36.5xt -yt，z' t = = xt yt- zt，x0 = = z0 = = 0，y0= =1，x， y，z，t，0，20， MaxSteps3000Out3=x InterpolatingFunction0. ，20.，< >，y In

15、terpolatingFunction0.，20.， < >，zInterpolatingFunction0.，20.， < >，ln4: =ParametricPlot3DEvaluatext，yt，zt / . %，t，0，20，PlotPoints 1000图13-45 3维相轨线Out3= -Graphics3D-说明：以上范例中ln1取高精度，而且是复系数方程。ln2是求解在x=1时的近似值,1i1i求精确解能得到准确值e4，读者可以求e4的近似值与Out2的结果比较，验证近似解的精确 度确实很高。ln3在求解时增大步数，成功地得到了由In 4绘制的如图13-

16、45所示的解的相 轨线。ln4所示的绘图语句与前面例子中的不同，现在只要会模仿使用它们就行了，要想弄 清原理请参阅相关Mathematica书籍。三、拉氏变换Mathematica可以进行拉普拉斯变换，其变换使用的函数调用格式如下：Lap laceTra nsformf，t， s求函数f (t)的Lap lace变换，返回自变量为 s的函数。InverseLaplaceTransformF, s，t求函数 F(s)的 Laplace逆变换，返回自变量为t的函数。其中函数f (t)和F (s)也可以是函数表，这样可一次变换多个函数 例8求函数t 4和et sint的拉氏变换。解：In1: =La

17、placeTransformL4，t，s24Out1=飞sln2: =LaplaceTransformExpt Sint，t，sOut2=12 - 2s s2ln3: =lnverseLaplaceTransform%1，s，t4Out3=t4ln4: =lnverseLaplaceTransform%2, s，tOut4= 1 iie(1_i)t(_1 - e2it)2ln5: =FullSimplify%Out 5=et Si nt例9 求函数f (t)= t3 eat的拉氏变换解：ln1: = LaplaceTransformQ3 Expa t，t，sOut1=6(-as)4以上只是直接

18、进行拉氏变换和逆变换的例子。以下使用拉氏变换解常微分方程，解法原 理见本书理论篇，这里完全实现了计算机求解。例10用拉氏变换解微分方程：x + 3x + 3x' + x = 1 满足条件 x (0) = x ' (0) = x(0) = 0 的解。解:ln1: =f仁LaplaceTransform xt+3x t+ 3x ' t+xt,t, sOut1=LaplaceTransformxt , t, s +3s LaplaceTransformxt, t, s +23 (sLaplaceTransformxt, t, s - x0) - s x0 +2 f3 (s L

19、aplaceTransformxt, t, s - s x0 - x ' 0) s x' 0- x 0ln2: =s仁LaplaceTransform1, t, s1s(1 s)Out2=ln3:=x 0= x ' 0= x0=0 ;Solvef1= =s1, LaplaceTransformxt, t, sOut4=LaplaceTra nsform xt, t, s >3In5:s, t1=I nv erseLaplaceTra nsform-s(1 + s)1 _L2Out5= 1 ex(2 2t t2)说明：上例中的LaplaceTransformxt,

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

21、 + y' t - 2xt,t, s;In2: =s1=LaplaceTransform0, - 2Exp-t , t, s;ln3: = x0=3 ; x' 0= 2 ; y0=0 ;Solvef1= =s1 , LaplaceTransformxt , t, s,LaplaceTransformyt, t, s;ln5: =lnverseLaplaceTransformFlatte nLaplaceTransformxt , t, s,LaplaceTransformyt, t, s / . % , s, tOut5=5 - e-t - 3et + 2e2t, e-t (-

22、 1 + et) 2 (1 + 2 et) In6: =Simplify%-t t2t -t t 2tOut6= 5 - e - 3e + 2e , e - 3e + 2e 说明：在上例中，不显示任何中间结果，语句比较简练。其中， ln1和ln2分别对方程 组的左边和右边进行拉氏变换，ln3解出X(s)和丫(s)0 In5比较难懂，可以参看前面 的例题，这里是从Out3中自动将解X (s)和丫(s)提取出来，再进行拉氏逆变换。Out5是x (t), y (t) , Out6将答案化简。本例已经将求解过程一般化，只需改变方程组和初 值的数据，就可以解其它方程组了。四、级数1.求和与求积求有限或无

23、穷和、积的函数是：i maxSumf, i , imin , imax求二 f (i)，其中 imin 可以是-x, minimax可以是x(即+x)，但是必须满足imin < imax。基本输入模板中也有求和专用的符号， 使用模板输入更方便。求多重和，也可以使用基本输入模i max求I丨f (i),基本输入模板中也有i=imin求多重积，也可以使用基本输入Sumf, i , imin , imax , j , jmin, jmax， 板连续多次输入求和符号得到。Productf , i , imin, imax求积符号Productf , i , imin , imax, j , jm

24、in , jmax,模板连续多次输入求积符号得到。例12求下列级数的和与积:n(1) k2 ,(2)k =1二 1J,(3)' 丄,(4)k=1 kO0 1II ek2。k=183解：In1: =SumkA2 , k , 1 , nOut1= 1 n(1n)(1 2n)Out2=6In2:oO='1你八2k =1712QOln3: = ' 1/ kk-1Sum： div : Sum does not conv erge.:1Out3=' -k 4 kcdln4: =j【Exp1/kA2kJOut4= e莎说明：上例中第三个级数发散，Mathematica给出提示

25、，并在不能给出结果时将输入的式 子作为输出。NSum和NProduct得到数值解。2. 将函数展开为幕级数将函数展开为幕级数的函数调用格式如下:Seriesf, x , xo, n为止。Seriesf, x , xo, n, y , yo, m例13展开下列函数为幕级数：sin x(1) y=tgx,( 2) y 匸x将函数f (x)在xo处展成幕级数直到n次项将函数f (x, y)先对y后对x展开。 y = f (x),( 4) y = eTyo解：ln1:=SeriesTanx,x ,0 , 935x2x17x762x9Out1=x3153152835ox10ln2: =SeriesSin

26、x /x , x , 0 , 92468亠x x xx10Out2= 1ox61205040362880ln3: =Seriesfx , x , 1 , 71(x-1)31 1OUt3= f1 )2f 1)(1)2 61241720f 1(x-1)6f 1(x-1)75 0 4 0ox -18ln4: =SeriesExpx y, x , 0 , 3 , y , 0 , 2/ 2、Out4= 1 (y oy3)x - iy oy3 x2 oy3x3 ox4l2丿说明：上例中ln3表明也可以展开抽象的函数。对已经展开的幕级数进行操作的两个函数是：Normalexpr将幕级数expr去掉余项转换成

27、多项式SeriesCoefficientexpr, n找出幕级数 expr 的 n 次项系数。例14将y = arcsinx展开为幕级数，只取前9项并去掉余项。解：In1: =SeriesArcSinx，x，0，9Out1=3579L x ± 3x ± 5x i 35x x640112115210oxIn2:=Normal%Out2=x33x55x7 35x9=x亠亠亠640112 1152In3: =SeriesCoefficient%1，53Out3= 403.傅里叶级数求傅里叶级数就是求出傅里叶系数，傅里叶系数是一个积分表达式，所以利用积分函数In tegrate就可以实现。例如，设周期矩形脉冲信号的脉冲宽度为t,脉冲幅度为E,周期