Mathematica求解方程组级数.docx

Mathematica求解方程组级数.docx

《Mathematica求解方程组级数.docx》由会员分享，可在线阅读，更多相关《Mathematica求解方程组级数.docx（12页珍藏版）》请在冰豆网上搜索。

Mathematica求解方程组级数

方程（组）与级数的Mathematica求解

[学习目标]

1.能用Mathematica求各种方程（组）的数值解和近似解；

2.能对常见函数进行幂级数的展开。

一、求解简单方程（组）

数学里的方程是带有变量的等式。

一般地说，一个或一组方程总是对于方程中出现的变量的可能取值范围增加了一些限制。

所谓求解方程就是设法把方程对于变量取值的限制弄清楚，最好的结果是用不含变量的表达式把变量的值表示出来。

在这个系统里，方程也用含有变量的等式表示，要注意的是在这里等号用连续的两个等号（==）表示。

方程的两端可以是任何数学表达式。

用户可以自己操作Mathematica系统去求解方程，例如使用移项一类的等价变换规则对方程加以变形、对方程的两端进行整理、把函数作用于方程的两端等等。

系统也提供了一些用于求解方程的函数。

1、求方程的代数解

最基本的方程求解函数是Solve，它可以用于求解方程（主要是多项式方程）或方程组。

Solve有两个参数，第一个参数是一个方程，或者是由若干个方程组的表（表示一个方程组）;第二个参数是要求解的变量或变量表。

例如，下面的式子对于变量X求解方程

：

In[1]:

=Solve[x^4-x^3-6x^2+1==0,x]

输入了这个表达式，系统立刻就能计算出方程的四个根，求出的解都是精确解（代数根）。

对于一般的多项式，这样得出的解常常是用根式描述的复数。

方程的解被表示成一个表，表中是几个子表，每一个子表的形式都是{x->...}，箭头后面是方程的一个解。

Solve也可以求解多变量的方程或者方程组:

In[2]:

=Solve[{x-2y==0,x^2-y==1},{x,y}]

这个表达式求解方程组:

有时求解方程会得到非常复杂的解。

例如将上面的第一个方程稍加变形，所得到的解的表达式就会变得很长:

In[3]:

=Solve[x^4-x^3-6x^2=2==0,x]

这个表达式求出的解的表达式非常长，以至一个计算机屏幕显示不下。

使用MS-DOS系统上的Mathematica的读者可以用键盘上的PgUP键和PgDn键把计算机屏幕上已经卷出的表达式翻回来阅读，附录B里提供了使用这类计算机的有关操作的更详细的说明。

对于使用图形界面提供的功能去翻阅前面的结果。

在被求解的方程里还可以有其他符号参数，可以要求系统对于这一个或者那一个变量求解方程。

对于Mathematica系统来说，方程中的符号变量（无论使用什么变量名）都是一样的。

对于处理复杂的方程，MATHEMETICA系统还提供了例外两个有用的函数。

函数Eliminate用于从方程组消去一个或几个变量，例如下面的表达式消去方程组里的变量Y:

IN[4]:

=Eliminate[{X^2-2Y==1,X+2Y==4},Y]

Eliminate的使用形式与Solve类似，它的第二个参数用于说明希望消去的变量。

另一个函数Reduce用于化简复杂的方程或方程组，它试图用一组比较简单的逻辑关系来描述由原来方程所描述的变量之间的关系。

它的使用形式与Solve,Eliminate一样，这里不举例字了。

2、求方程的数值解

理论上已经证明，对于五次以上的多项式方程没有求代数解的一般方法，MATHEMATICA也求不出那些不能分解因式的五次以上的多项式方程的解，例如：

IN[5]:

=SOLVE[X^5+5X^3-2==0,X]

它返回一个带有函数TORUOES的表达式。

可以把函数N作用到这个结果表达式上，求出方程的数值解：

IN[6]:

=N[%]

可以看到系统同时求出了方程的五个根的时候可以直接用函数N和SOLVE结合完成工作：

IN[7]:

=N[SOLVE[X^6+4X^2-31==0,X]]

在系统里直接提供了一个函数NSOLVE做这件事。

对于更复杂的方程（或方程组），用SOLVE求不出根，使用函数N也解决不了问题。

对于这样的方程，用户可以使用REDUCE,ELIMINATE等函数去处理，设法把方程描述的变量之间的关系搞清楚。

如果需要的就是方程的根，那么只要用求数值根的函数FINDROOT。

函数FINDROOT求数值根所采用的方法与人们一般用计算机求数值根的方法一样。

但是，由于MATHEMATIC有求导函数的能力，在这里计算有导函数的表达式的数值根就非常简单。

不管表达式多么复杂，系统都能自动的求出它的导函数。

求数值根使用的也是牛顿法，用户必须给FINDROOT提供一个初始值。

下面一个简单的例子：

IN[8]:

=FINDROOT[SIN[X]EXP[2X]-COS[X]==0,{X,0.5}]

对于求不出导函数的表达式，例如用户自己定义的一个复杂计算函数，使用FindRoot提供函数值取不同符号（正负号）的两个点（用表的形式放在上面初始值0.5的位置），形式是：

IN[9]:

=FindRoot[FUN1[X]==0,{X,{0,1}}]

这里假使FUN1是用户定义的一个函数。

使用计算机求数值根的第一个问题是确定初始点，若初始值选取得不好将给求根带来困难。

再一个麻烦是用户要自己求出函数的导函数。

在MATHEMATICA系统里处理的对象是表达式，一个表达式可以服务于不同的用途，可以作为求值的对象，作为画图的对象，也可以作为演算的对象。

当需要求一个表达式的数值根的时候，表达式的这样的多种功

能，或者说MATHEMATICA系统对于表达式的多方面的操作能力就表现出很大的优越性。

一个代数表达式，无论多么复杂，MATHEMATICA系统都可以直接求出它的导函数，可以作出它的图形。

从图形上我们很容易认识这个函数表达式在某一个区间的大致性质，包括它的根的出现和分布情况。

对表达式的这些认识为人们确定如何取初始值、如何求根提供了很有价值的线索。

这样，某些比较难以处理的问题可能就容易解决了。

二、求解常微分方程（组）

1、常微分方程（组）的精确解

Mathematica能求常微分方程（组）的准确解，能求解的类型大致覆盖了人工求解的范围，功能很强。

但不如人灵活（例如在隐函数和隐方程的处理方面），输出的结果与教材上的答案可能在形式上不同。

另外，Mathematica求数值解也很方便，且有利于作出解的图形。

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

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求微分方程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]→

}}

2、常微分方程（组）的数值解

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

NDSolve[eqns，{y1，y2，…}，{x，xmin，xmax}]求常微分方程（组）的近似解。

其中微分方程和初值条件的表示法如同DSolve，未知函数仍有带自变量和不带自变量两种形式，通常使用后一种更方便。

初值点x0可以取在区间[xmin，xmax]上的任何一点处，得到插值函数InterpolatingFunction[domain，table]类型的近似解，近似解的定义域domain一般为[domain，table]，也有可能缩小。

例4求常微分方程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所示。

例5求常微分方程组：

满足初始条件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。

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

（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书籍。

三、级数

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}，…]求多重积，也可以使用基本输入模板连续多次输入求积符号得到。

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展开。

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

（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次项系数。

例8将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]=

习题

1．求下列方程（组）的通解或特解。

（1）x^3-5x+1==0，

（2）x^6-5x^4+x^3-2==0

（3）

（4）

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

（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．求下列级数的和或积：

（1）

，

（2）

，（3）

，（4）

4.将下列函数展开为x的幂级数。

（1）

；

（2）y=cos2x；

（3）y=（1-x）ln（1-x）