1、第四章 方程求解教学目的：学习并掌握计算机代数系统Maple下进行代数方程和微分方程求解的方法和技巧，并了解其在方程求解中的缺限。 教学目标：掌握代数方程和微分方程求解的方法和技巧，并尝试应用所学数学基础解决Maple下关于方程求解的缺陷问题。 重点内容：代数方程求解，微分方程求解。 难点内容：高次代数方程求解，非线性微分方程求解。1 代数方程(组)求解1.1 常用求解工具solve 求解代数方程或代数方程组, 使用Maple中的solve函数. 求解关于 x的方程eqn=0的命令格式为： solve(eqn, x); 求解关于变量组vars的方程组eqns的命令为： solve(eqns,
2、vars); eqn:=(x2+x+2)*(x-1); solve(eqn,x); 当然, solve也可以求解含有未知参数的方程： eqn:=2*x2-5*a*x=1; solve(eqn,x); solve函数的第一个参数是有待求解的方程或方程的集合, 当然也可以是单个表达式或者表达式的集合, 如下例： solve(a+ln(x-3)-ln(x),x);对于第二个参数, Maple的标准形式是未知变量或者变量集合 , 当其被省略时, 函数indets 自动获取未知变量. 但当方程中含有参数时, 则会出现一些意想不到的情况： solve(a+ln(x-3)-ln(x);很多情况下, 我们知道
3、一类方程或方程组有解, 但却没有解决这类方程的一般解法, 或者说没有解析解.比如,一般的五次或五次以上的多项式,其解不能写成解析表达式.Maple具备用所有一般算法尝试所遇到的问题,在找不到解的时候, Maple 会用RootOf 给出形式解. x7-2*x6-4*x5-x3+x2+6*x+4; solve(%); solve(cos(x)=x,x);但是, 有时候, Maple甚至对一些“ 显而易见”的结果置之不理, 如： solve(sin(x)=3*x/Pi,x); 此方程的解为0 ,6, 但Maple 却对这个超越方程无能为力 , 即便使用allvalues求解也只有下述结果： all
4、values(%); 另外一个问题是, Maple在求解方程之前,会对所有的方程或表达式进行化简, 而不管表达式的类型, 由此而产生一些低级的错误： (x-1)2/(x2-1); solve(%);1 但是, 大量实验表明 , solve的确是一个实用的方程求解工具, 但是也不可盲目相信它给出的一切结果, 特别是对于非线性方程而言, 对于给出的结果需要加以验证. 下面通过几个例子说明在Maple中非线性方程组的求解问题. 例：求解方程组： eqns:=x2+y2=25,y=x2-5; vars:=x,y; solve(eqns,vars);也可用下面的语句一步求出： solve(x2+y2=2
5、5,y=x2-5,x,y); 这个问题非常简单, 但通常遇到的非线性问题却不是这么简单, 例如要求解方程组： eqns:=x2+y2=1,sqrt(x+y)=x-y; vars:=x,y; sols:=solve(eqns,vars); 可以看出, 方程解的形式是以集合的序列给出的, 序列中的每一个集合是方程的一组解, 这样就很利于我们用subs把解代入原方程组进行检验： subs(sols2,eqns); sols2:=allvalues(sols1); simplify(subs(sols2,eqns);1.2 其他求解工具1.2.1 数值求解对于求代数方程的数值解问题, Maple提供了
6、函数fsolve, fsolve的使用方法和solve 很相似: fsolve(eqns, vars, options); 其中, eqns表示一个方程、方程组或者一个程序, vars表示一个未知量或者未知量集合, options控制解的参数(诸如：complex : 复根; maxsols=n：只找到n阶最小根; intervals：在给定闭区间内求根, 等)。 fsolve(x5-x+1,x);1.167303978 fsolve(x5-x+1,x,complex); fsolve(x3-3*x+1,x,0.1);.3472963553 对于多项式方程, fsolve在默认情况下以给出所有
7、的实数解,如果附加参数complex, 就可以给出所有的解。但对于更一般的其他形式的方程, fsolve却往往只满足于得到一个解： eqn:=sin(x)=x/2; fsolve(eqn);0. fsolve(eqn,x,0.1.infinity); 1.895494267 fsolve(eqn,x,-infinity.-0.1); 1.895494267函数fsolve主要基于两个算法, 通常使用牛顿法，如果牛顿法无效, 它就改而使用切线法。为了使fsolve可以求得所有的实根，我们通常需要确定这些根所在的区间. 对于单变量多项式，函数realroot 可以获得多项式的所有实根所在的区间.
8、4+6*x+x2-x3-4*x5-2*x6+x7; realroot(%);函数realroot还有一个可选参数 , 它是用来限制区间的最大长度的 , 为了保证使用数值求解方法时收敛, 我们可以用它限制区间的最大长度： realroot(%,1/1000);求解方程或方程组的整数解时使用函数isolve, 它常常被用来求解不定方程. 例如著名的“ 百钱买百鸡”问题的求解过程为： isolve(x+y+z=100,5*x+3*y+z/3=100);当_Z1=1,2,3时，据此可得满足该问题的三组解为： x, y, z=4, 18, 78, x, y, z=8, 11, 81, x, y, z=1
9、2, 4, 84 1.2.2 整数环中的方程(组)求解利用Maple中的函数msolve(eqns, vars, n), 可以在模n的整数环中求解方程(组)eqns . 例：在Z 7中求解Pell方程 7328yx msolve(y7=x3-28,7); 再如下例： msolve(y4=x3+32,5); 1.2.3 递归方程的求解在Maple中, 可以求解有限差分方程 (也称递归方程 ), 所需调用的函数是rsolve, 该函数使用的是一些比较通用的方法, 例如产生函数法、z变换法以及一些基于变量替换和特征方程的方法. 作为例子, 求解Fibonacci多项式: eq:=f(n)=f(n-1
10、)+2*f(n-2); rsolve(eq,f(0)=1,f(1)=1,f(n); 当然, 并不是所有的递归形式的函数方程的解可以写成解析形式, 如果不能, Maple将保留原来的调用形式。此时，可用asympt函数获得它的渐进表达式, 也就是1/n的级数解。例如，对于一个具有超越形式的递归函数方程, 仍然可以得到解的渐进形式： rsolve(u(n+1)=ln(u(n)+1),u(n); asympt(%,n,5); 1.2.4 不等式(组)求解求解一元不等式方程(组)使用命令solve： solve(x-1)*(x-2)*(x-3) solve(x-1+a)*(x-2+a)*(x-3+a)
11、 solve(exp(x)x+1); solve(x2*y2=0,x-y=1,x with(simplex): cnsts:=3*x+4*y-3*z obj:=-x+y+2*z; maximize(obj,cnsts union x=0,y=0,z=0); 但是，客观地讲，Maple中求解不等式(组)是一个薄弱环节，尤其对非线性不等式组几乎无能为力。2 常微分方程求解微分方程求解是数学研究与应用的一个重点和难点。Maple能够显式或隐式地解析地求解许多微分方程求解。在常微分方程求解器dsolve中使用了一些传统的技术例如laplace变换和积分因子法等，函数pdesolve 则使用诸如特征根法
12、等经典方法求解偏微分方程。此外，Maple还提供了可作摄动解的所有工具，例如Poincare-Lindstedt法和高阶多重尺度法。帮助处理常微分方程(组)的各类函数存于Detools软件包中，函数种类主要有：可视化类的函数，处理宠加莱动态系统的函数，调整微分方程的函数，处理积分因子、李对称法和常微分方程分类的函数，微分算子的函数，利用可积性与微分消去的方法简化微分方程的函数，以及构造封闭解的函数等。更重要的是其提供的强大的图形绘制命令Deplot能够帮助我们解决一些较为复杂的问题。2.1 常微分方程的解析解求解常微分方程最简单的方法是利用求解函数dsolve. 命令格式为： dsolve(O
13、DE); dsolve(ODE, y(x), extra_args); dsolve(ODE, ICs, y(x), extra_args); dsolve(sysODE, ICs, funcs,extra_args); 其中，ODE 常微分方程，y(x)单变量的任意变量函数，Ics初始条件，sysODEODE方程组的集合，funcs 变量函数的集合，extra_args 依赖于要求解的问题类型。 例如，对于一阶常微分方程 (可用dsolve直接求得解析解： lnxyy ODE:=x*diff(y(x),x)=y(x)*ln(x*y(x)-y(x); dsolve(ODE,y(x); 可以看出
14、，dsolve的第一个参数是待求的微分方程，第二个参数是未知函数。需要注意的是，无论在方程中还是作为第二个参数，未知函数必须用函数的形式给出(即：必须加括号, 并在其中明确自变量 )，这一规定是必须的，否则Maple将无法区分方程中的函数、自变量和参变量，这一点和我们平时的书写习惯不一致。为了使其与我们的习惯一致，可用alias 将函数用别称表示： alias(y=y(x); ODE:=x*diff(y,x)=y*ln(x*y)-y; dsolve(ODE,y); 函数dsolve给出的是微分方程的通解, 其中的任意常数是用下划线起始的内部变量表示的. 在Maple中, 微分方程的解是很容易验
15、证的 , 只需要将解代入到原方程并化简就可以了. subs(%,ODE); assume(x,real): assume(_C1,real):simplify(%); evalb(%);evalb函数的目的是对一个包含关系型运算符的表达式使用三值逻辑系统求值，返回的值是true, false和FAIL。如果无法求值，则返回一个未求值的表达式。通常包含关系型运算符“=, =”的表达式在Maple中看作是代数方程或者不等式。然而，作为参数传递给evalb或者出现在 if或while语句的逻辑表达式中时，它们会被求值为true或false 。值得注意的是，evalb不化简表达式，因此在使用evalb
16、之前应将表达式化简，否则可能会出错。再看下面常微分方程的求解：12+= yy alias(y=y(x): ODE:=diff(y,x)=sqrt(y2+1); dsolve(ODE,y); 函数dsolve对于求解含有未知参变量的常微分方程也完全可以胜任: alias(y=y(x): ODE:=diff(y,x)=-y/sqrt(a2-y2); sol:=dsolve(ODE,y); 由此可见，对于不能表示成显式结果的微分方程解，Maple 尽可能将结果表示成隐式解。另外，对于平凡解 y=0 常常忽略，这一点应该引起注意。dsolve对于求解微分方程初值问题也十分方便的: ODE:=diff(
17、u(t),t$2)+omega2*u(t)=0; dsolve(ODE,u(0)=u0,D(u)(0)=v0,u(t); 2.2 利用积分变换求解微分方程对于特殊的微分方程，我们还可以指定dsolve利用积分变换方法求解，只需要在dsolve中加入可选参数method=transform即可。其中transform是积分变换，可以是laplace 、fourier、fouriercos或者fouriersin 变换。作为例子, 我们来看一个具有阻尼的振子在阶跃冲击(Heaviside函数)下的响应： ODE:=diff(u(t),t$2)+2*d*omega*diff(u(t),t)+omeg
18、a2*u(t)=Heaviside(t); initvals:=(u(0)=u0,D(u)(0)=v0); solution:=dsolve(ODE,initvals,u(t),method=laplace); Maple给出了问题的通解, 但没有区分自由振动(d=0 )、欠阻尼(01 )的情况. 下面加以区分求解： assume(omega0): simplify(subs(d=0,solution); K:=subs(d=1/5,u0=1,v0=1,solution); with(plots):plot3d(rhs(%),omega=2/3.4/3,t=0.20,style=hidden,
19、orientation=-30,45,axes=framed);对于d=1的情况，可用下式获得结果： limit(rhs(solution),d=1);再如下例： diff(u(t),t$2)+3*diff(u(t),t)+2*u(t)=exp(-abs(t); dsolve(%,u(t),method=fourier); 2.3 常微分方程组的求解函数dsolve不仅可以用来求解单个常微分方程，也可以求解联立的常微分方程组. 特别是对于线性微分方程组，由于数学上具有成熟的理论，Maple 的求解也是得心应手。其命令格式为： dsolve( eqn1, eqn2, , ini_conds, v
20、ars); 其中, ini_conds 是初始条件. eqn1:=diff(x(t),t)=x(t)+y(t),diff(y(t),t)=y(t)-x(t); dsolve(eqn1,x(t),y(t); eqn2:=2*diff(x(t),t$2)+2*x(t)+y(t)=2*t; eqn3:=diff(y(t),t$2)+2*x(t)+y(t)=t2+1; dsolve(eqn2, qn3, x(0)=0, D(x)(0)=1, y(0)=0, D(y)(0)=0, x(t),y(t) ); = ()xt + + 18()sin2t2112t3148t434t, = ()yt + + 14
21、()sin2t212t12t216t3124t4 2.4 常微分方程的级数解法1) 泰勒级数解法当一个常微分方程的解析解难以求得时，可以用Maple求得方程解的级数近似，这在大多数情况下是一种非常好的方法。级数解法是一种半解析半数值的方法。泰勒级数法的使用命令为： dsolve(ODE,Ics, y(x), series); 或dsolve(ODE,Ics, y(x), type=series); 下面求解物理摆的大幅振动方程： l , 其中l是摆长, 是摆角, g是重力加速度. ODE:=l*diff(theta(t),t$2)=-g*sin(theta(t); initvals:=theta(0)=0,D(theta)(0)=v0/l; sol:=dsolve(ODE,initvals,theta(t),type=series); Order:=11: sol:=dsolve(ODE,initvals,theta(t),type=series);
