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

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

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

完整版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．将函数