完整版Mathematica常微分方程拉氏变换与级数实验.docx

完整版Mathematica常微分方程拉氏变换与级数实验.docx

《完整版Mathematica常微分方程拉氏变换与级数实验.docx》由会员分享，可在线阅读，更多相关《完整版Mathematica常微分方程拉氏变换与级数实验.docx（16页珍藏版）》请在冰豆网上搜索。

完整版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解下列常微分方程（组）：

（1）

，

（2）

，（3）

，

（4）

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

解：

In[1]：

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

y[x]，x]

Out[1]=

In[2]：

=DSolve[y′[x]==（1+y[x]^2）/（（x+x^3）y[x]），y[x]，x]

Out[2]={{

}，{

}}

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）

，

（2）

，（3）

。

解：

In[1]：

=DSolve[

+3y″[x]+3y′[x]+y[x]==（x-5）Exp[-x]，

y[x]，x]

Out[1]={{

}}

In[2]：

=Simplify[%]

Out[2]={{

}}

In[3]：

=DSolve[x^2+y′[x]^2==1，y[x]，x]

Out[3]={{

}，

{

}}

In[4]：

=DSolve[Sqrt[y′[x]]==xy[x]，y[x]，x]

Out[4]={{

}}

说明：

由以上可以看出对方程的类型并无限制，但是输出的答案未必符合习惯，例如第一个方程的答案需要化简，有时即使化简后也未必与教材上的答案一致。

例3求微分方程

的通解。

解：

In[1]：

=DSolve[y′[x]+2xy[x]==xE^（-x^2），y[x]，x]

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

}}

这就是所给微分方程的通解。

式中的C[1]是通解中的任意常数。

上述命令也可以输入为：

DSolve[D[y[x]]+2xy[x]==xE^（-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]

Out[1]={{y[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]表明返回的解放在一个表中，不便使用，实际的解就是插值函数：

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

In[2]的结果是用y表示解函数的名字，因此In[3]顺利画出解曲线如图13-43所示。

例6求常微分方程组：

满足初始条件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最大步数。

StartingStepSize初始步长。

以上可选参数的默认值都为Automatic，其中AccuracyGoal和PrecisionGoal的默认值比WorkingPrecision小10，当解趋于0时应将AccuracyGoal取成Infinity。

对于常微分方程，最大步长默认值为1000。

这个函数也可以解偏微分方程，最大步长默认值为200。

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

（1）

，满足初始条件y（0）=1的特解；

（2）

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

解：

In[1]：

=NDSolve[{y′[x]==I/4y[x]，y[0]==1}，y，{x，1}，

AccuracyGoal→20，PrecisionGoal→20，WorkingPrecision→25]

Out[1]={{y→InterpolatingFunction[

{{0，1.000000000000000000000000000}}，<>]}

In[2]：

=y[1]/.%

Out[2]={0.968912424710644784118519+

0.2474039592545229296234109ⅱ}

In[3]：

=NDSolve[{x′[t]==-3（x[t]-y[t]），

y′[t]==-x[t]z[t]+36.5x[t]-y[t]，

z′[t]==x[t]y[t]-z[t]，

x[0]==z[0]==0，y[0]==1}，

{x，y，z}，{t，0，20}，MaxSteps→3000]

Out[3]={{x→InterpolatingFunction[{{0.，20.}}，<>]，

y→InterpolatingFunction[{{0.，20.}}，<>]，

z→InterpolatingFunction[{{0.，20.}}，<>]}}，

In[4]：

=ParametricPlot3D[Evaluate[{x[t]，y[t]，z[t]}/.%]，

{t，0，20}，PlotPoints→1000]

图13-453维相轨线

Out[3]=-Graphics3D-

说明：

以上范例中In[1]取高精度，而且是复系数方程。

In[2]是求解在x=1时的近似值，求精确解能得到准确值

，读者可以求

的近似值与Out[2]的结果比较，验证近似解的精确度确实很高。

In[3]在求解时增大步数，成功地得到了由In[4]绘制的如图13-45所示的解的相轨线。

In[4]所示的绘图语句与前面例子中的不同，现在只要会模仿使用它们就行了，要想弄清原理请参阅相关Mathematica书籍。

三、拉氏变换

Mathematica可以进行拉普拉斯变换，其变换使用的函数调用格式如下：

LaplaceTransform[f，t，s]求函数f（t）的Laplace变换，返回自变量为s的函数。

InverseLaplaceTransform[F，s，t]求函数F（s）的Laplace逆变换，返回自变量为t的函数。

其中函数f（t）和F（s）也可以是函数表，这样可一次变换多个函数。

例8求函数t4和etsint的拉氏变换。

解：

In[1]：

=LaplaceTransform[t^4，t，s]

Out[1]=

In[2]：

=LaplaceTransform[Exp[t]Sin[t]，t，s]

Out[2]=

In[3]：

=InverseLaplaceTransform[%1，s，t]

Out[3]=t4

In[4]：

=InverseLaplaceTransform[%2，s，t]

Out[4]=

In[5]：

=FullSimplify[%]

Out[5]=etSin[t]

例9求函数f（t）=t3eat的拉氏变换。

解：

In[1]：

=LaplaceTransform[t^3Exp[at]，t，s]

Out[1]=

以上只是直接进行拉氏变换和逆变换的例子。

以下使用拉氏变换解常微分方程，解法原理见本书理论篇，这里完全实现了计算机求解。

例10用拉氏变换解微分方程：

+3x″+3x′+x=1满足条件x″（0）=x′（0）=x（0）=0的解。

解：

In[1]：

=f1=LaplaceTransform[

[t]+3x″[t]+3x′[t]+x[t]，

t，s]

Out[1]=LaplaceTransform[x[t]，t，s]+

s3LaplaceTransform[x[t]，t，s]+

3（sLaplaceTransform[x[t]，t，s]-x[0]）-s2x[0]+

3（s2LaplaceTransform[x[t]，t，s]-sx[0]-x′[0]）-

sx′[0]-x″[0]

In[2]：

=s1=LaplaceTransform[1，t，s]

Out[2]=

In[3]：

=x″[0]=x′[0]=x[0]=0；

Solve[f1==s1，LaplaceTransform[x[t]，t，s]]

Out[4]={{

}}

In[5]：

=InverseLaplaceTransform[

，s，t]

Out[5]=

说明：

上例中的LaplaceTransform[x[t]，t，s]就是教材中的X（s），In[3]解出X（s），其余过程与教科书完全相同。

现在可以将一切计算留给计算机，学生只要弄清解法原理及过程。

技巧：

充分利用复制、粘贴功能，可以加快输入速度，避免键入错误。

上例中In[5]就可以从Out[4]中将表达式复制过来。

例11求微分方程组：

满足条件x（0）=3，x′（0）=2，y（0）=0的特解。

解：

In[1]：

=f1=LaplaceTransform[

{x″[t]-2x′[t]-y′[t]+2y[t]，x′[t]+y′[t]-2x[t]}，

t，s]；

In[2]：

=s1=LaplaceTransform[{0，-2Exp[-t]}，t，s]；

In[3]：

=x[0]=3；x′[0]=2；y[0]=0；

Solve[{f1==s1}，{LaplaceTransform[x[t]，t，s]，

LaplaceTransform[y[t]，t，s]}]；

In[5]：

=InverseLaplaceTransform[

Flatten[

{LaplaceTransform[x[t]，t，s]，

LaplaceTransform[y[t]，t，s]}/.%]，s，t]

Out[5]={5-e-t-3et+2e2t，e-t（-1+et）2（1+2et）}

In[6]：

=Simplify[%]

Out[6]={5-e-t-3et+2e2t，e-t-3et+2e2t}

说明：

在上例中，不显示任何中间结果，语句比较简练。

其中，In[1]和In[2]分别对方程组的左边和右边进行拉氏变换，In[3]解出X（s）和Y（s）。

In[5]比较难懂，可以参看前面的例题，这里是从Out[3]中自动将解X（s）和Y（s）提取出来，再进行拉氏逆变换。

Out[5]是{x（t），y（t）}，Out[6]将答案化简。

本例已经将求解过程一般化，只需改变方程组和初值的数据，就可以解其它方程组了。

四、级数

1．求和与求积

求有限或无穷和、积的函数是：

Sum[f，{i，imin，imax}]求

，其中imin可以是-∞，imax可以是∞（即+∞），但是必须满足imin≤imax。

基本输入模板中也有求和专用的符号，使用模板输入更方便。

Sum[f，{i，imin，imax}，{j，jmin，jmax}，…]求多重和，也可以使用基本输入模板连续多次输入求和符号得到。

Product[f，{i，imin，imax}]求

，基本输入模板中也有求积符号。

Product[f，{i，imin，imax}，{j，jmin，jmax}，…]求多重积，也可以使用基本输入模板连续多次输入求积符号得到。

例12求下列级数的和与积：

（1）

，

（2）

，（3）

，（4）

。

解：

In[1]：

=Sum[k^2，{k，1，n}]

Out[1]=

In[2]：

=

Out[2]=

In[3]：

=

Sum：

：

div：

Sumdoesnotconverge.

Out[3]=

In[4]：

=

Out[4]=

说明：

上例中第三个级数发散，Mathematica给出提示，并在不能给出结果时将输入的式子作为输出。

NSum和NProduct得到数值解。

2．将函数展开为幂级数

将函数展开为幂级数的函数调用格式如下：

Series[f，{x，x0，n}]将函数f（x）在x0处展成幂级数直到n次项为止。

Series[f，{x，x0，n}，{y，y0，m}]将函数f（x，y）先对y后对x展开。

例13展开下列函数为幂级数：

（1）y=tgx，

（2）

，（3）y=f（x），（4）y=exy。

解：

In[1]：

=Series[Tan[x]，{x，0，9}]

Out[1]=

In[2]：

=Series[Sin[x]/x，{x，0，9}]

Out[2]=

In[3]：

=Series[f[x]，{x，1，7}]

Out[3]=

In[4]：

=Series[Exp[xy]，{x，0，3}，{y，0，2}]

Out[4]=

说明：

上例中In[3]表明也可以展开抽象的函数。

对已经展开的幂级数进行操作的两个函数是：

Normal[expr]将幂级数expr去掉余项转换成多项式。

SeriesCoefficient[expr，n]找出幂级数expr的n次项系数。

例14将y=arcsinx展开为幂级数，只取前9项并去掉余项。

解：

In[1]：

=Series[ArcSin[x]，{x，0，9}]

Out[1]=

In[2]：

=Normal[%]

Out[2]=

In[3]：

=SeriesCoefficient[%1，5]

Out[3]=

3．傅里叶级数

求傅里叶级数就是求出傅里叶系数，傅里叶系数是一个积分表达式，所以利用积分函数Integrate就可以实现。

例如，设周期矩形脉冲信号的脉冲宽度为τ，脉冲幅度为E，周期为T，这种信号在一个周期[

，

]内的表达式为

求其傅里叶级数时，可以先求出傅利叶系数。

为了和Mathematica中的常数E相区分，以下用Ee表示脉冲幅度，用tao表示脉冲宽度τ，根据傅利叶系数的积分表达式，输入以下语句：

a0=2/TIntegrate[Ee，

{t，-tao/2，tao/2}]

a[n_]=2/TIntegrate[EeCos[2nPit/T]，

{t，-tao/2，tao/2}]

b[n_]=2/TIntegrate[EeSin[2nPit/T]，

{t，-tao/2，tao/2}]

可得到下面三个输出，即分别是a0，an与bn，即

a0=

，an=

与bn=0

从而可写出给定的傅利叶级数为：

习题13.5

1．求下列一阶微分方程的通解或特解。

（1）y′-3xy=2x；

（2）xy′+y-ex=0，y|x=a=b。

2．求下列二阶微分方程的通解或特解。

（1）y″-2y′+5y=5x+2；

（2）y″+2y′+2y=-e-x，y（0）=y′（0）=0。

3．用拉氏变换解微分方程组：

4．展开下列函数为x幂级数，并求其收敛区间。

；

（2）y=cos2x；（3）y=（1-x）ln（1-x）

5．将函数

分别展开成正弦级数和余弦级数。

6．求下列函数的拉氏变换：

（1）t2+6t-3；

（2）e-4tsin3tcos2t；（3）

。

7．求下列像函数的拉氏逆变换：

（1）F（p）=

；

（2）F（p）=

；

（3）F（p）=

。