发布时间 : 星期一 文章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[x0]= =y0}，y[x]，x] 求满足初始条件y（x0）= y0的特解y（x）。
DSolve[{eqn1，eqn2，…}，{y1[x]，y2[x]，…}，x] 求方程组的通解。 DSolve[{equ1，…，y1[x0]= =y10，…}，{y1[x]，y2[x]，…}，x] 求方程组的特解。 说明：应当特别注意，方程及各项参数的表述方式很严格，容易出现输入错误。微分方程的表示法只有通过例题才能说清楚。
例1 解下列常微分方程（组）：
2y1?y22?(x?1)，（2）y?? （1）y??， （3） x?1(x?x3)y5常微分方程、拉氏变换与级数实验
?y??z， ??z??y??y??z （4）?的通解及满足初始条件y（0）=0，z（0）=1的特解。
?z???y 解：In[1]：=DSolve[y′[x]= =2y[x]/（x+1）+（x+1）^（5/2）， y[x]，x]
??2?? Out[1]=??y[x]?(1?x)7/2?(1?x)2c[1]??
3???? In[2]：=DSolve[y′[x]= =（1+y[x]^2）/((x+x^3）y[x])，y[x]，x]
?1? Out[2]={{y[x]??11?c[1]?1??c[1]22xx}， {y[x]?}} 111?21?2xx 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）y????3y???3y??y?(x?5)e?x，（2）x2?(y?)2?1，（3）y??xy。 解：In[1]：=DSolve[y???[x]+3y″[x] +3y′[x] + y[x] = =（x - 5）Exp[-x]， y[x]，x]
1?x2?x2??x?5x2x3? Out[1]={{y[x]?ex??5x???ex???? ????22?3???25x3x4??x?x?x2 1e?x??????eC[1]?exC[2]?exC[3]}} ??2?34? In[2]：=Simplify[%] Out[2]={{y[x]?1?xe(?20x3?x4?24C[1]?24xC[2]?24x2C[3])}} 24 In[3]：=DSolve[x^2 + y′[x]^2 = = 1，y[x]，x]
1ArcSin[x]?C[1]}， Out[3]={{y[x]??x1?x2?22 {y[x]?1ArcSin[x]x1?x2??C[1]}} 22 In[4]：=DSolve[Sqrt[y′[x]] = = x y[x]，y[x]，x] Out[4]={{y[x]??3}} 3x?C[1] 说明：由以上可以看出对方程的类型并无限制，但是输出的答案未必符合习惯，例如第一个方程的答案需要化简，有时即使化简后也未必与教材上的答案一致。
例3 求微分方程
2dy?2xy?xe?x的通解。 dx 解：In[1]：=DSolve[y′[x]+2x y[x]= = x E^（-x^2），y[x]，x]
66
212 Out[1]={{y[x]→e?xx2?e?xC[1]}}
2 这就是所给微分方程的通解。式中的C[1]是通解中的任意常数。
上述命令也可以输入为：DSolve[D[y[x]] + 2x y[x]= =x E^（ - x^2），y[x]，x]。
例4 求微分方程xy′+ y - ex = 0在初始条件y|x=1 = 2e下的特解。 解：In[1]：=DSolve[x*y′[x]+y[x]-E^x= =0，y[1]= =2E，y[x]，x]
e?ex
Out[1]= {{y[x]→}}
x
二、 常微分方程（组）的数值解
函数NDSolve用于求给定初值条件或边界条件的常微分方程（组）的近似解，其调用格式如下：
NDSolve[eqns，{y1，y2，…}，{x，xmin，xmax}] 求常微分方程（组）的近似解。 其中微分方程和初值条件的表示法如同DSolve，未知函数仍有带自变量和不带自变量两种形式，通常使用后一种更方便。初值点x0可以取在区间[xmin，xmax]上的任何一点处，得到插值函数InterpolatingFunction[domain，table]类型的近似解，近似解的定义域domain一般为[domain，table]，也有可能缩小。
例5 求常微分方程y′= x2 + y2，满足初始条件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]表明返回的解放在一个表中，不便使用，实际的解就是插值函数：
67
InterpolatingFunction[{{-2.，2.}}，< >]
In[2]的结果是用y表示解函数的名字，因此In[3]顺利画出解曲线如图13-43所示。
例6 求常微分方程组：
1???x?y?x3?x ? 3?y???x? 满足初始条件x（0）=0，y（0）=1的数值解。 解：In[1]：=s1=NDSolve[{x′[t]= = y[t] -（x[t]^3/3 - x[t]）， y′[t]= = - x[t]，x[0]= =0，y[0]= =1}， {x，y}，{t，-15，15}]
Out[1]={{x→InterpolatingFunction[{{-15.，15.}}，< >]， y→InterpolatingFunction[{{-15.，15.}}，< >]}} In[2]：= x=x / . s1[[1，1]] y=y / . s1[[1，2]]
Out[2]=InterpolatingFunction[{{-15.，15.}}，< >] Out[3]=InterpolatingFunction[{{-15.，15.}}，< >] In[4]：=ParametricPlot[{x[t]，y[t]}，{t，-15，15}， AspectRatio→Automatic]
图13-44 解的相轨线
Out[3]= -Graphics-
说明：上例是求一个著名方程组的近似解，其中In[2]也可以改用一个赋值式{x，y}={x，y} / . Flatten[s1]，一次得到两个函数。通过求数值解容易得到它的相图，In[4]绘制了解的相轨线如图13-44所示，图中表明原点是奇点，极限环的形状也已经得到。 为了应付复杂的情况，需要设置可选参数：
WorkingPrecision 参见数值积分部分的介绍。 AccuracyGoal 计算结果的绝对误差。 PrecisionGoal 计算结果的相对误差。 MaxSteps 最大步数。
68