(完整版)Mathematica——常微分方程、拉氏变换与级数实验

§13.5

常微分方程、拉氏变换与级数实验 [学习目标]

1. 会用Mathematica 求解微分方程（组）；

2. 能用Mathematica 求微分方程（组）的数值解；

3. 会利用Mathematica 进行拉氏变换与逆变换；

4. 能进行幂级数和傅里叶级数的展开。

一、 常微分方程（组）

Mathematica 能求常微分方程（组）的准确解，能求解的类型大致覆盖了人工求解的范围，功能很强。但不如人灵活（例如在隐函数和隐方程的处理方面），输出的结果与教材上的答案可能在形式上不同。另外，Mathematica 求数值解也很方便，且有利于作出解的图形。在本节中，使用Laplace 变换解常微分方程（组）的例子也是十分成功的，过去敬而远之的方法如今可以轻而易举的实现了。

求准确解的函数调用格式如下：

DSolve[eqn ，y[x]，x] 求方程eqn 的通解y （x ），其中自变量是x 。

DSolve[{eqn ，y[x 0]= =y 0}，y[x]，x] 求满足初始条件y （x 0）= y 0的特解y （x ）。

DSolve[{eqn1，eqn2，…}，{y 1[x]，y 2[x]，…}，x] 求方程组的通解。 DSolve[{equ1，…，y 1[x 0]= =y 10，…}，{y 1[x]，y 2[x]，…}，x] 求方程组的特解。 说明：应当特别注意，方程及各项参数的表述方式很严格，容易出现输入错误。微分方程的表示法只有通过例题才能说清楚。

例1 解下列常微分方程（组）：

（1）25)1(12+++='x x y y ，（2）y x x y y )(132++='， （3） ⎩

⎨⎧-='='y z z y ， （4）⎩⎨⎧-='='y

z z y 的通解及满足初始条件y （0）=0，z （0）=1的特解。

解：In[1]：=DSolve[y ′[x]= =2y[x]/（x+1）+（x+1）^（5/2），

y[x]，x]

Out[1]=⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩

⎨⎧+++→]1[)1()1(32][22/7c x x x y In[2]：=DSolve[y ′[x]= =（1+y[x]^2）/((x+x^3）y[x])，y[x]，x]

Out[2]={{2211]1[11][x c x x y ++-

--→}， {2211]1[11][x

c x x y ++--→}}

In[3]：=DSolve[{y ′[x]= =z[x]，z ′[x]= = -y[x]}，

{y[x]，z[x]}，x]

Out[3]={{y[x]→C[1]Cos[x]+ C[2]Sin[x]，

z[x]→C[2]Cos[x]- C[1]Sin[x]}}

In[4]：=DSolve[{y ′[x]= =z[x]，z ′[x]= = -y[x]，y[0]= =0，z[0]= =1}，

{y[x]，z[x]}，x]

Out[4]={{y[x]→Sin[x]，z[x]→Cos[x]}}

提示：认真观察上例，可以从中学习输入格式，未知函数总带有自变量，等号用连续键入两个等号表示，这两点由于不习惯会出错！导数符号用键盘上的撇号，连续两撇表示二阶导数，这与习惯相同。自变量、未知量、初始值的表示法与普通变量相同。

说明：输出结果总是尽量用显式解表出，有时反而会使表达式变得复杂，这与教科书的习惯不同。当求显式解遇到问题时，会给出提示。通解中的任意常数用C[1]，C[2]，…表示。

例2 求解下列微分方程：

（1）x e x y y y y --=+'+''+''')5(33，（2）1)(22='+y x ，（3）xy y ='。

解：In[1]：=DSolve[][x y '''+3y ″[x] +3y ′[x] + y[x] = =（x - 5）Exp[-x]，

y[x]，x]

Out[1]={{+⎪⎪⎭⎫ ⎝⎛-+⎪⎪⎭⎫ ⎝⎛+-→--32

52521][3222x x x e x x x e x y x x ]3[]2[]1[43521243C x e xC e C e x x e x x x x ----+++⎪⎪⎭

⎫ ⎝⎛+-}} In[2]：=Simplify[%]

Out[2]={{])3[24]2[24]1[2420(24

1][243C x xC C x x e x y x ++++-→-}} In[3]：=DSolve[x^2 + y ′[x]^2 = = 1，y[x]，x]

Out[3]={{]1[2

][121][2C x ArcSin x x x y +---→}， {]1[2

][121][2C x ArcSin x x x y +--→}} In[4]：=DSolve[Sqrt[y ′[x]] = = x y[x]，y[x]，x]

Out[4]={{]

1[3][3C x x y --→}} 说明：由以上可以看出对方程的类型并无限制，但是输出的答案未必符合习惯，例如第一个方程的答案需要化简，有时即使化简后也未必与教材上的答案一致。

例3 求微分方程22x xe xy dx

dy -=+的通解。 解：In[1]：=DSolve[y ′[x]+2x y[x]= = x E^（-x^2），y[x]，x]

Out[1]={{y[x]→]1[2

1222C e x e x x --+}} 这就是所给微分方程的通解。式中的C[1]是通解中的任意常数。

上述命令也可以输入为：DSolve[D[y[x]] + 2x y[x]= =x E^（ - x^2），y[x]，x]。

例4 求微分方程xy ′+ y - e x = 0在初始条件y|x=1 = 2e 下的特解。

解：In[1]：=DSolve[x*y ′[x]+y[x]-E^x= =0，y[1]= =2E ，y[x]，x]

Out[1]= {{y[x]→x

e e x

+}} 二、 常微分方程（组）的数值解

函数NDSolve 用于求给定初值条件或边界条件的常微分方程（组）的近似解，其调用格式如下：

NDSolve[eqns ，{y 1，y 2，…}，{x ，xmin ，xmax}] 求常微分方程（组）的近似解。 其中微分方程和初值条件的表示法如同DSolve ，未知函数仍有带自变量和不带自变量两种形式，通常使用后一种更方便。初值点x 0可以取在区间[xmin ，xmax]上的任何一点处，得到插值函数InterpolatingFunction[domain ，table]类型的近似解，近似解的定义域domain 一般为[domain ，table]，也有可能缩小。

例5 求常微分方程y ′= x 2 + y 2，满足初始条件y （0）= 0的数值解。

解：In[1]：=s1=NDSolve[{y ′[x]= =x^2+y[x]^2，y[0]= =0}，

y ，{x ，-2，2}]

Out[1]={{y →InterpolatingFunction[{{-2.，2.}}，< >]}}

In[2]：= y=y / . s1[[1]]

Out[2]=InterpolatingFunction[{{-2.，2.}}，< >]

In[3]：=Plot[y[x]，{x ，-2，2}，AspectRatio →Automatic ，

PlotRange →{-1.5，1.5}]

图13-43 微分方程的解曲线

Out[3]= -Graphics-

上例中包含许多值得学习的实用内容，其中第二项参数使用y 而不是y[x]，比用y[x]好。如果求解区间改为{x ，-3，3}，就会出现警告提示，实际得不到[-3，3]上的解。Out[1]表明返回的解放在一个表中，不便使用，实际的解就是插值函数：