Mathematica——常微分方程、拉氏变换与级数实验.doc

13.5 常微分方程、拉氏变换与级数实验学习目标1. 会用Mathematica求解微分方程（组）；2. 能用Mathematica求微分方程（组）的数值解；3. 会利用Mathematica进行拉氏变换与逆变换；4. 能进行幂级数和傅里叶级数的展开。一、 常微分方程（组） Mathematica能求常微分方程（组）的准确解，能求解的类型大致覆盖了人工求解的范围，功能很强。但不如人灵活（例如在隐函数和隐方程的处理方面），输出的结果与教材上的答案可能在形式上不同。另外，Mathematica求数值解也很方便，且有利于作出解的图形。在本节中，使用Laplace变换解常微分方程（组）的例子也是十分成功的，过去敬而远之的方法如今可以轻而易举的实现了。 求准确解的函数调用格式如下： DSolveeqn，yx，x 求方程eqn的通解y（x），其中自变量是x。 DSolveeqn，yx0= =y0，yx，x 求满足初始条件y（x0）= y0的特解y（x）。 DSolveeqn1，eqn2，y1x，y2x，x 求方程组的通解。 DSolveequ1，y1x0= =y10，y1x，y2x，x 求方程组的特解。 说明：应当特别注意，方程及各项参数的表述方式很严格，容易出现输入错误。微分方程的表示法只有通过例题才能说清楚。例1 解下列常微分方程（组）： （1），（2）， （3） ， （4）的通解及满足初始条件y（0）=0，z（0）=1的特解。 解：In1：=DSolveyx= =2yx/（x+1）+（x+1）（5/2）， yx，x Out1= In2：=DSolveyx= =（1+yx2）/(x+x3）yx)，yx，x Out2=， In3：=DSolveyx= =zx，zx= = -yx， yx，zx，x Out3=yxC1Cosx+ C2Sinx， zxC2Cosx- C1Sinx In4：=DSolveyx= =zx，zx= = -yx，y0= =0，z0= =1， yx，zx，x Out4=yxSinx，zxCosx 提示：认真观察上例，可以从中学习输入格式，未知函数总带有自变量，等号用连续键入两个等号表示，这两点由于不习惯会出错！导数符号用键盘上的撇号，连续两撇表示二阶导数，这与习惯相同。自变量、未知量、初始值的表示法与普通变量相同。 说明：输出结果总是尽量用显式解表出，有时反而会使表达式变得复杂，这与教科书的习惯不同。当求显式解遇到问题时，会给出提示。通解中的任意常数用C1，C2，表示。例2 求解下列微分方程： （1），（2），（3）。解：In1：=DSolve+3yx +3yx + yx = =（x - 5）Exp-x， yx，x Out1= In2：=Simplify% Out2= In3：=DSolvex2 + yx2 = = 1，yx，x Out3=， In4：=DSolveSqrtyx = = x yx，yx，x Out4= 说明：由以上可以看出对方程的类型并无限制，但是输出的答案未必符合习惯，例如第一个方程的答案需要化简，有时即使化简后也未必与教材上的答案一致。例3 求微分方程的通解。 解：In1：=DSolveyx+2x yx= = x E（-x2），yx，x Out1=yx 这就是所给微分方程的通解。式中的C1是通解中的任意常数。 上述命令也可以输入为：DSolveDyx + 2x yx= =x E（ - x2），yx，x。例4 求微分方程xy+ y - ex = 0在初始条件y|x=1 = 2e下的特解。 解：In1：=DSolvex*yx+yx-Ex= =0，y1= =2E，yx，x Out1= yx二、 常微分方程（组）的数值解 函数NDSolve用于求给定初值条件或边界条件的常微分方程（组）的近似解，其调用格式如下： NDSolveeqns，y1，y2，x，xmin，xmax 求常微分方程（组）的近似解。 其中微分方程和初值条件的表示法如同DSolve，未知函数仍有带自变量和不带自变量两种形式，通常使用后一种更方便。初值点x0可以取在区间xmin，xmax上的任何一点处，得到插值函数InterpolatingFunctiondomain，table类型的近似解，近似解的定义域domain一般为domain，table，也有可能缩小。例5 求常微分方程y= x2 + y2，满足初始条件y（0）= 0的数值解。解：In1：=s1=NDSolveyx= =x2+yx2，y0= =0， y，x，-2，2 Out1=yInterpolatingFunction-2.，2.， In2：= y=y / . s11 Out2=InterpolatingFunction-2.，2.， In3：=Plotyx，x，-2，2，AspectRatioAutomatic， PlotRange-1.5，1.5图13-43 微分方程的解曲线 Out3= -Graphics- 上例中包含许多值得学习的实用内容，其中第二项参数使用y而不是yx，比用yx好。如果求解区间改为x，-3，3，就会出现警告提示，实际得不到-3，3上的解。Out1表明返回的解放在一个表中，不便使用，实际的解就是插值函数： InterpolatingFunction-2.，2.， In2的结果是用y表示解函数的名字，因此In3顺利画出解曲线如图13-43所示。例6 求常微分方程组： 满足初始条件x（0）=0，y（0）=1的数值解。 解：In1：=s1=NDSolvext= = yt -（xt3/3 - xt）， yt= = - xt，x0= =0，y0= =1， x，y，t，-15，15 Out1=xInterpolatingFunction-15.，15.， yInterpolatingFunction-15.，15.， In2：= x=x / . s11，1 y=y / . s11，2 Out2=InterpolatingFunction-15.，15.， Out3=InterpolatingFunction-15.，15.， In4：=ParametricPlotxt，yt，t，-15，15， AspectRatioAutomatic图13-44 解的相轨线 Out3= -Graphics- 说明：上例是求一个著名方程组的近似解，其中In2也可以改用一个赋值式x，y=x，y / . Flattens1，一次得到两个函数。通过求数值解容易得到它的相图，In4绘制了解的相轨线如图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的特解。 解：In1：=NDSolveyx= =I/4yx，y0= =1，y，x，1， AccuracyGoal20，PrecisionGoal20，WorkingPrecision25 Out1=yInterpolatingFunction 0，1.000000000000000000000000000， In2：=y1 / . % Out2=0.968912424710644784118519 + 0.2474039592545229296234109 In3：=NDSolvext= = -3（xt -yt）， yt = = -xt zt+36.5xt -yt， zt = = xt yt- zt， x0 = = z0 = = 0，y0= =1， x，y，z，t，0，20，MaxSteps3000 Out3=xInterpolatingFunction0.，20.， yInterpolatingFunction0.，20.， zInterpolatingFunction0.，20.， In4：=ParametricPlot3DEvaluatext，yt，zt / . %， t，0，20，PlotPoints 1000图13-45 3维相轨线 Out3= -Graphics3D-说明：以上范例中In1取高精度，而且是复系数方程。In2是求解在x=1时的近似值，求精确解能得到准确值，读者可以求的近似值与Out2的结果比较，验证近似解的精确度确实很高。In3在求解时增大步数，成功地得到了由In4绘制的如图13-45所示的解的相轨线。In4所示的绘图语句与前面例子中的不同，现在只要会模仿使用它们就行了，要想弄清原理请参阅相关Mathematica书籍。三、 拉氏变换 Mathematica可以进行拉普拉斯变换，其变换使用的函数调用格式如下： LaplaceTransformf，t，s 求函数f（t）的Laplace变换，返回自变量为s的函数。 InverseLaplaceTransformF，s，t 求函数F（s）的Laplace逆变换，返回自变量为t的函数。 其中函数f（t）和F（s）也可以是函数表，这样可一次变换多个函数。例8 求函数t 4和et sint的拉氏变换。 解：In1：=LaplaceTransformt4，t，s Out1= In2：=LaplaceTransformExpt Sint，t，s Out2= In3：=InverseLaplaceTransform%1，s，t Out3=t4 In4：=InverseLaplaceTransform%2，s，t Out4= In5：=FullSimplify% Out5=et Sint例9 求函数f（t）= t3 eat的拉氏变换。 解：In1：= LaplaceTransformt3 Expa t，t，s Out1= 以上只是直接进行拉氏变换和逆变换的例子。以下使用拉氏变换解常微分方程，解法原理见本书理论篇，这里完全实现了计算机求解。例10 用拉氏变换解微分方程： + 3x+ 3x+ x = 1满足条件x(0) = x(0) = x(0) = 0的解。解：In1：=f1=LaplaceTransform t +3xt+ 3xt+xt， t，s Out1=LaplaceTransformxt，t，s + s3LaplaceTransformxt，t，s + 3（sLaplaceTransformxt，t，s - x0）- s2x0 + 3（s2LaplaceTransformxt，t，s - s x0 - x0）- s x0- x0 In2：=s1=LaplaceTransform1，t，s Out2= In3：= x0= x0= x0=0； Solvef1= =s1，LaplaceTransformxt，t，s Out4= In5：=InverseLaplaceTransform，s，t Out5= 说明：上例中的LaplaceTransformxt，t，s就是教材中的X（s），In3解出X（s），其余过程与教科书完全相同。现在可以将一切计算留给计算机，学生只要弄清解法原理及过程。 技巧：充分利用复制、粘贴功能，可以加快输入速度，避免键入错误。上例中In5就可以从Out4中将表达式复制过来。例11 求微分方程组： 满足条件x（0）=3，x（0）=2，y（0）=0的特解。解：In1：=f1=LaplaceTransform xt - 2xt - yt + 2yt，xt + yt - 2xt， t，s； In2：=s1=LaplaceTransform0，- 2Exp-t，t，s； In3：= x0=3； x0= 2；y0=0； Solvef1= =s1，LaplaceTransformxt，t，s， LaplaceTransformyt，t，s； In5：=InverseLaplaceTransform Flatten LaplaceTransformxt，t，s， LaplaceTransformyt，t，s / . %，s，t Out5=5 - e-t - 3et + 2e2t，e-t（- 1 + et）2（1 + 2 et） In6：=Simplify% Out6= 5 - e-t - 3et + 2e2t，e-t - 3et + 2e2t 说明：在上例中，不显示任何中间结果，语句比较简练。其中，In1和In2分别对方程组的左边和右边进行拉氏变换，In3解出X（s）和Y（s）。In5比较难懂，可以参看前面的例题，这里是从Out3中自动将解X（s）和Y（s）提取出来，再进行拉氏逆变换。Out5是x（t），y（t），Out6将答案化简。本例已经将求解过程一般化，只需改变方程组和初值的数据，就可以解其它方程组了。四、 级数1 求和与求积 求有限或无穷和、积的函数是： Sumf，i，imin，imax 求，其中imin可以是-，imax可以是（即+），但是必须满足iminimax。基本输入模板中也有求和专用的符号，使用模板输入更方便。 Sumf，i，imin，imax，j，jmin，jmax， 求多重和，也可以使用基本输入模板连续多次输入求和符号得到。 Productf，i，imin，imax 求，基本输入模板中也有求积符号。 Productf，i，imin，imax，j，jmin，jmax， 求多重积 ，也可以使用基本输入模板连续多次输入求积符号得到。例12 求下列级数的和与积： （1），（2） ，（3） ，（4） 。 解：In1：=Sumk2，k，1，n Out1= In2：= Out2= In3：= Sum：div：Sum does not converge. Out3= In4：= Out4= 说明：上例中第三个级数发散，Mathematica给出提示，并在不能给出结果时将输入的式子作为输出。 NSum和NProduct得到数值解。2 将函数展开为幂级数将函数展开为幂级数的函数调用格式如下： Seriesf，x，x0，n 将函数f（x）在x0 处展成幂级数直到n次项为止。 Seriesf，x，x0，n，y，y0，m 将函数f（x，y）先对y后对x展开。例13 展开下列函数为幂级数：（1） y=tgx，（2） ， （3）y = f（x），（4）y = exy。 解：In1：=SeriesTanx，x，0，9 Out1= In2：=SeriesSinx /x，x，0，9 Out2= In3：=Seriesfx，x，1，7 Out3= In4：=SeriesExpx y，x，0，3，y，0，2 Out4= 说明：上例中In3表明也可以展开抽象的函数。 对已经展开的幂级数进行操作的两个函数是： Normalexpr 将幂级数expr去掉余项转换成多项式。 SeriesCoefficientexpr，n 找出幂级数expr的n次项系数。例14 将y = arcsinx展开为幂级数，只取前9项并去掉余项。 解：In1：=SeriesArcSinx，x，0，9 Out1= In2：=Normal% Out2= In3：=SeriesCoefficient