计算机代数系统第4章-方程求解

第四章 方程求解

1 代数方程(组)求解

1.1 常用求解工具—solve

求解代数方程或代数方程组, 使用Maple中的solve函数. 求解关于x的方程eqn=0的命令格式为：

solve(eqn, x);

求解关于变量组vars的方程组eqns的命令为： solve(eqns, vars);

> eqn:=(x^2+x+2)*(x-1);

eqn := (x???x???2)(x???1)2

> solve(eqn,x);

1,?12???12I7,?12???12I7

当然, solve也可以求解含有未知参数的方程： > eqn:=2*x^2-5*a*x=1;

eqn := 2x???5ax???12

> solve(eqn,x);

51a???4425a???8,254a???1425a???8

2solve函数的第一个参数是有待求解的方程或方程的集合, 当然也可以是单个表达

式或者表达式的集合, 如下例：

> solve(a+ln(x-3)-ln(x),x);

3eaa?1???e

对于第二个参数, Maple的标准形式是未知变量或者变量集合, 当其被省略时, 函数indets自动获取未知变量. 但当方程中含有参数时, 则会出现一些意想不到的情况： > solve(a+ln(x-3)-ln(x));

- 100 -

{x???x,a????ln(x???3)???ln(x)}

很多情况下, 我们知道一类方程或方程组有解, 但却没有解决这类方程的一般解法, 或者说没有解析解. 比如, 一般的五次或五次以上的多项式, 其解不能写成解析表达式. Maple具备用所有一般算法尝试所遇到的问题, 在找不到解的时候, Maple会用RootOf给出形式解.

> x^7-2*x^6-4*x^5-x^3+x^2+6*x+4;

x???2x???4x???x???x???6x???476532

> solve(%);

1???5,1???5,RootOf(_Z5???_Z???1,index???1),RootOf(_Z5???_Z???1,index???2),RootOf(_Z5???_Z???1,index???3),RootOf(_Z5???_Z???1,index???4),RootOf(_Z5???_Z???1,index???5)

> solve(cos(x)=x,x);

RootOf(_Z???cos(_Z))

对于方程组解的个数可用nops命令获得, 如： > eqns:={seq(x[i]^2=x[i],i=1..7)};

eqns := {x1???x1,x2???x2,x3???x3,x4???x4,x5???x5,x6???x6,x7???x7}

2222222> nops({solve(eqns)});

128

但是, 有时候, Maple甚至对一些“显而易见”的结果置之不理, 如： > solve(sin(x)=3*x/Pi,x);

RootOf(3_Z???sin(_Z)?)

此方程的解为??6, 0, 但Maple却对这个超越方程无能为力, 即便使用allvalues

求解也只有下述结果： > allvalues(%);

RootOf(3_Z???sin(_Z)?,0.)

另外一个问题是, Maple在求解方程之前,会对所有的方程或表达式进行化简, 而不管表达式的类型, 由此而产生一些低级的错误： > (x-1)^2/(x^2-1);

(x???1)x???122

> solve(%);

1

- 101 -

但是, 大量实验表明, solve的确是一个实用的方程求解工具, 但是也不可盲目相信它给出的一切结果, 特别是对于非线性方程而言, 对于给出的结果需要加以验证.

下面通过几个例子说明在Maple中非线性方程组的求解问题. ?x2?y2?25例：求解方程组：?2

?x?9?y> eqns:={x^2+y^2=25,y=x^2-5};

eqns := {y???x???5,x???y???25}222

> vars:={x,y};

vars := {x,y}

> solve(eqns,vars);

{x???0,y???-5},{x???0,y???-5},{y???4,x???3},{y???4,x???-3}

也可用下面的语句一步求出：

> solve({x^2+y^2=25,y=x^2-5},{x,y});

{x???0,y???-5},{x???0,y???-5},{y???4,x???3},{y???4,x???-3}

这个问题非常简单, 但通常遇到的非线性问题却不是这么简单, 例如要求解方程组：x?y?1,22x?y?x?y

> eqns:={x^2+y^2=1,sqrt(x+y)=x-y}; vars:={x,y};

eqns := {x???y???1,22x???y???x???y}

vars := {x,y}

> sols:=solve(eqns,vars);

sols := {y???RootOf(2_Z2???4_Z???3,-1.000000000x????RootOf(2_Z2???4_Z???3,-1.000000000???.7071067812???.7071067812I),I)???2},{x???1,y???0}

可以看出, 方程解的形式是以集合的序列给出的, 序列中的每一个集合是方程的一组解, 这样就很利于我们用subs把解代入原方程组进行检验： > subs(sols[2],eqns);

{1???1}

> sols2:=allvalues(sols[1]);

11sols2 := {x????1???I2,y????1???I22> simplify(subs(sols2,eqns));

- 102 -

2}

{I2???I2,1???1}

1.2 其他求解工具

1.2.1 数值求解

对于求代数方程的数值解问题, Maple提供了函数fsolve, fsolve的使用方法和solve很相似:

fsolve(eqns, vars, options);

其中, eqns表示一个方程、方程组或者一个程序, vars表示一个未知量或者未知量集合, options控制解的参数(诸如：complex: 复根; maxsols=n：只找到n阶最小根; intervals：在给定闭区间内求根, 等). > fsolve(x^5-x+1,x);

-1.167303978

I,.7648844336???.3524715460I,> fsolve(x^5-x+1,x,complex);

-1.167303978,-.1812324445???1.083954101II,-.1812324445???1.083954101.7648844336???.3524715460

> fsolve(x^3-3*x+1,x,0..1);

.3472963553

对于多项式方程, fsolve在默认情况下以给出所有的实数解, 如果附加参数complex, 就可以给出所有的解. 但对于更一般的其他形式的方程, fsolve却往往只满足于得到一个解：

> eqn:=sin(x)=x/2;

1eqn := 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可以获得多项式的所有实根所在的区间. > 4+6*x+x^2-x^3-4*x^5-2*x^6+x^7;

4???6x???x???x???4x???2x???x23567

> realroot(%);

- 103 -

[[0,2],[2,4],[-2,-1]]

函数realroot还有一个可选参数, 它是用来限制区间的最大长度的, 为了保证使用数值求解方法时收敛, 我们可以用它限制区间的最大长度： > realroot(%%,1/1000);

??1195,299?,?3313,1657?,?-633,-1265??????1024256????1024512????5121024???? ????????求解方程或方程组的整数解时使用函数isolve, 它常常被用来求解不定方程. 例如

著名的“百钱买百鸡”问题?的求解过程为：

> isolve({x+y+z=100,5*x+3*y+z/3=100});

{z???75???3_Z1,x???4_Z1,y???25???7_Z1}

据此可得满足该问题的三组解为：

{x, y, z}={4, 18, 78}, {x, y, z}={8, 11, 81}, {x, y, z}={12, 4, 84}

1.2.2 整数环中的方程(组)求解

利用Maple中的函数msolve(eqns, vars, n), 可以在模n的整数环中求解方程(组)eqns.

例：在Z7中求解Pell方程y7?x3?28 > msolve(y^7=x^3-28,7);

{x???3,y???6},{x???4,y???1},{y???0,x???0},{x???1,y???1},{y???6,x???6},{x???2,y???1},{y???6,x???5}

再如下例：

> msolve(y^4=x^3+32,5);

{x???2,y???0},{x???4,y???1},{x???4,y???2},{x???4,y???3},{x???4,y???4}

1.2.3 递归方程的求解

在Maple中, 可以求解有限差分方程(也称递归方程), 所需调用的函数是rsolve, 该函数使用的是一些比较通用的方法, 例如产生函数法、z变换法以及一些基于变量替换和特征方程的方法. 作为例子, 求解Fibonacci多项式: > eq:=f(n)=f(n-1)+2*f(n-2);

eq := f(n)???f(n???1)???2f(n???2)

> rsolve({eq,f(0)=1,f(1)=1},f(n));

12nn(-1)???2 33当然, 并不是所有的递归形式的函数方程的解可以写成解析形式, 如果不能, Maple将保留原来的调用形式. 此时, 可用asympt函数获得它的渐进表达式, 也就是1/n的级

数解. 例如, 对于一个具有超越形式的递归函数方程, 仍然可以得到解的渐进形式： ?

百钱买百鸡问题：用100元钱买100只鸡, 大公鸡5元钱1只, 大母鸡3元钱1只, 小鸡1元钱3只, 问如何买法？

- 104 -

> rsolve(u(n+1)=ln(u(n)+1),u(n));

rsolve(u(n???1)???ln(u(n)???1),u(n))

> asympt(%,n,5);

21n_C??????23n2ln(n)???19???13_C???122?222_C?????_C???ln(n)???ln(n)????29?91??3???O??4?3?n?n??

1.2.4 不等式(组)求解

求解一元不等式方程(组)使用命令solve： > solve((x-1)*(x-2)*(x-3)<0,x);

RealRange(??,Open(1)),RealRange(Open(2),Open(3))

> solve((x-1+a)*(x-2+a)*(x-3+a) < 0, {x});

{x???1???a},{2???a???x,x???3???a}

> solve(exp(x)>x+1);

RealRange(??,Open(0)),RealRange(Open(0),?)

> solve({x^2*y^2=0,x-y=1,x<>0});

{y???0,x???1},{y???0,x???1}

对于由不等式方程组约束的最优问题的求解使用“线性规则”工具包simplex:

> with(simplex):

> cnsts:={3*x+4*y-3*z<=23, 5*x-4*y-3*z<=10,7*x+4*y+11*z<=30};

cnsts := {3x???4y???3z???23,5x???4y???3z???10,7x???4y???11z???30}

> obj:=-x+y+2*z;

obj := ?x???y???2z

> maximize(obj,cnsts union {x>=0,y>=0,z>=0});

149{z???,y???,x???0}

282 常微分方程求解

微分方程求解是数学研究与应用的一个重点和难点. Maple能够显式或隐式地解析地求解许多微分方程求解. 在常微分方程求解器dsolve中使用了一些传统的技术例如laplace变换和积分因子法等, 函数pdesolve则使用诸如特征根法等经典方法求解偏微分

方程. 此外, Maple还提供了可作摄动解的所有工具, 例如Poincare-Lindstedt法和高阶多重尺度法.

帮助处理常微分方程(组)的各类函数存于Detools软件包中, 函数种类主要有：可视化类的函数, 处理宠加莱动态系统的函数, 调整微分方程的函数, 处理积分因子、李对称

- 105 -

法和常微分方程分类的函数, 微分算子的函数, 利用可积性与微分消去的方法简化微分方程的函数, 以及构造封闭解的函数等. 更重要的是其提供的强大的图形绘制命令Deplot能够帮助我们解决一些较为复杂的问题.

2.1 常微分方程的解析解

求解常微分方程最简单的方法是利用求解函数dsolve. 命令格式为： dsolve(ODE);

dsolve(ODE, y(x), extra_args);

dsolve({ODE, ICs}, y(x), extra_args);

dsolve({sysODE, ICs}, {funcs}, extra_args);

其中, ODE—常微分方程, y(x)—单变量的任意变量函数, Ics—初始条件, {sysODE}—ODE方程组的集合, {funcs}—变量函数的集合, extra_args—依赖于要求解的问题类型.

例如, 对于一阶常微分方程xy??yln(xy)?y可用dsolve直接求得解析解： > ODE:=x*diff(y(x),x)=y(x)*ln(x*y(x))-y(x);

?????y(x)ln(xy(x))???y(x)ODE := x? ???xy(x)????> dsolve(ODE,y(x));

e?x???_C1????y(x)???x

可以看出, dsolve的第一个参数是待求的微分方程, 第二个参数是未知函数. 需要

注意的是, 无论在方程中还是作为第二个参数, 未知函数必须用函数的形式给出(即：必须加括号, 并在其中明确自变量), 这一规定是必须的, 否则Maple将无法区分方程中的函数、自变量和参变量, 这一点和我们平时的书写习惯不一致. 为了使其与我们的习惯一致, 可用alias将函数用别称表示： > alias(y=y(x));

> ODE:=x*diff(y,x)=y*ln(x*y)-y;

??ODE := x????xy?????yln(xy)???y

??> dsolve(ODE,y);

y???e?x???_C1????x

函数dsolve给出的是微分方程的通解, 其中的任意常数是用下划线起始的内部变量

表示的.

在Maple中, 微分方程的解是很容易验证的, 只需要将解代入到原方程并化简就可以了.

- 106 -

> subs(%,ODE);

????_C1????ex???xx?x????????_C1???e??????x????????_C1???ln?exx????????_C1???e????xx????

> assume(x,real): assume(_C1,real): > simplify(%);

?e?x~???_C1~?????x~???_C1~????(?x~???_C1~)x~_C1~????e(?x~???_C1~)x~_C1~

> evalb(%);

true

evalb函数的目的是对一个包含关系型运算符的表达式使用三值逻辑系统求值, 返

回的值是true, false和FAIL. 如果无法求值, 则返回一个未求值的表达式. 通常包含关系型运算符“=, <>, <, <=, >, >=”的表达式在Maple中看作是代数方程或者不等式. 然而, 作为参数传递给evalb或者出现在if或while语句的逻辑表达式中时, 它们会被求值为true或false. 值得注意的是, evalb不化简表达式, 因此在使用evalb之前应将表达式化简, 否则可能会出错.

再看下面常微分方程的求解：y??y?1

2> alias(y=y(x)):

> ODE:=diff(y,x)=sqrt(y^2+1);

?2ODE := y???y???1

?x> dsolve(ODE,y);

y???sinh(x???_C1)

函数dsolve对于求解含有未知参变量的常微分方程也完全可以胜任: > alias(y=y(x)):

> ODE:=diff(y,x)=-y/sqrt(a^2-y^2);

ODE := ??xy????ya???y22

> sol:=dsolve(ODE,y);

?2a2???2aln??22?a???y???2aya22a???y22sol := x??????????_C1???0

由此可见, 对于不能表示成显式结果的微分方程解, Maple尽可能将结果表示成隐

- 107 -

式解. 另外, 对于平凡解y=0常常忽略, 这一点应该引起注意.

dsolve对于求解微分方程初值问题也十分方便的: > ODE:=diff(u(t),t$2)+omega^2*u(t)=0;

???2ODE := ?2u(t)?????u(t)???0??t???2

> dsolve({ODE,u(0)=u0,D(u)(0)=v0},u(t)); v0sin(?t)u(t)??????u0cos(?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)+omega^2*u(t)=Heaviside(t);

????2ODE := ?2u(t)????2d??u(t)?????u(t)???Heaviside??????t???t???2(t)

> initvals:=(u(0)=u[0],D(u)(0)=v[0]);

initvals := u(0)???u0,D(u)(0)???v0

> solution:=dsolve({ODE,initvals},u(t),method=laplace);

1?solution := u(t)??????e(?td?)?(?2u0???1)cosh(t?????d2?2????2)????(?v0???d?2u0???d)sinh(td2?2????2d2?2????2)?????

Maple给出了问题的通解, 但没有区分自由振动(d=0)、欠阻尼(01)的情况. 下面加以区分求解： > assume(omega>0):

> simplify(subs(d=0,solution));

1???cos(t?)?u(t)???2u0???cos(t?)???v0sin(t?)??2

> K:=subs(d=1/5,u[0]=1,v[0]=1,solution);

1?K := u(t)??????e(?1/5t?)???(?2???1)cosh?t????????????24??2??25?????????1?2???1?sinh??t????55????242??25?24??2??25????????????

> with(plots):

- 108 -

> plot3d(rhs(%%),omega=2/3..4/3,t=0..20,style=hidden,orientation=[-30,45],axes=framed);

对于d=1的情况, 可可用下式获得结果： > limit(rhs(solution),d=1);

(?u0????v0t???1????u0t???t????e?2223(t?))e(?t?)

再如下例：

> diff(u(t),t$2)+3*diff(u(t),t)+2*u(t)=exp(-abs(t));

????????3u(t)??t2???2??u(t)????2u(t)???e(?????t???t)

> dsolve(%,u(t),method=fourier);

u(t)???23e(?2t)Heaviside(t)???16eHeaviside(?t)???et(?t)tHeaviside(t)???12e(?t)Heaviside(t)

2.3 常微分方程组的求解

函数dsolve不仅可以用来求解单个常微分方程, 也可以求解联立的常微分方程组.

特别是对于线性微分方程组, 由于数学上具有成熟的理论, Maple的求解也是得心应手. 其命令格式为：

dsolve( {eqn1, eqn2, ?, ini_conds}, {vars}); 其中, ini_conds是初始条件.

> eqn1:={diff(x(t),t)=x(t)+y(t),diff(y(t),t)=y(t)-x(t)};

??eqn1 := {x(t)???x(t)???y(t),y(t)???y(t)???x(t)}

?t?t> dsolve(eqn1,{x(t),y(t)});

{x(t)???e(_C1sin(t)???_C2cos(t)),y(t)???e(_C1cos(t)???_C2sin(t))}tt

> eqn2:=2*diff(x(t),t$2)+2*x(t)+y(t)=2*t;

???eqn2 := 2?2x(t)????2x(t)???y(t)???2t??t???2

> eqn3:=diff(y(t),t$2)+2*x(t)+y(t)=t^2+1;

- 109 -

???2eqn3 := ?2y(t)????2x(t)???y(t)???t???1??t???2

> dsolve({eqn2, qn3, x(0)=0, D(x)(0)=1, y(0)=0, D(y)(0)=0}, {x(t),y(t)} );

{x(t)???18sin(2t)142???112t???12t???314812t???t???2434t,3y(t)???sin(2t)2???16t???124t}4

2.4 常微分方程的级数解法

1) 泰勒级数解法

当一个常微分方程的解析解难以求得时, 可以用Maple求得方程解的级数近似, 这

在大多数情况下是一种非常好的方法. 级数解法是一种半解析半数值的方法. 泰勒级数法的使用命令为：

dsolve({ODE,Ics}, y(x), 'series'); 或dsolve({ODE,Ics}, y(x), 'type=series');

下面求解物理摆的大幅振动方程：l???gsin?, 其中l是摆长, ?是摆角, g是重力加速度.

> ODE:=l*diff(theta(t),t$2)=-g*sin(theta(t));

???ODE := l?2?(t)?????gsin(?(t))??t???2

> initvals:=theta(0)=0,D(theta)(0)=v[0]/l;

initvals := ?(0)???0,D(?)(0)???v0l

> sol:=dsolve({ODE,initvals},theta(t),type=series);

sol := ?(t)???v0lt???1gv06l2t???31120gv0(v0???gl)l42t???O(t)56

> Order:=11:

> sol:=dsolve({ODE,initvals},theta(t),type=series);

sol := ?(t)???v0lt???1gv06l2t3???1120gv0(v0???gl)l42t5???15040gv0(11glv0???g2l2???v0)l624t7???1362880gv0(57gv0l???102g2v0l2???g3l3???v0)l8426t9???O(t11)

2) 幂级数解法

对于一个符号代数系统来说, 幂级数是必不可少的微分方程求解工具. 幂级数求解函数powsolve存于工具包powseries中. 但是, 这一求解函数的使用范围很有限, 它只可以用来求解多项式系数的线性常微分方程或方程组,其求解命令为：powseries[function] (prep)或直接载入软件包后用function(prep), prep为求解的线性微分方程及其初值.

- 110 -

例：求解：xy??y???4x2y?0

> ODE:=x*diff(y(x),x$2)+diff(y(x),x)+4*x^2*y(x)=0;

?????????4x2y(x)???0?????ODE := x?y(x)y(x)????x2???x????2

> dsolve(ODE,y(x));

4(3/2)??0,4x(3/2)?y(x)???_C1BesselJ?0,x???_C2BesselY?????3??3?

????> initvals:=y(0)=y0,D(y)(0)=0;

initvals := y(0)???y0,D(y)(0)???0

> with(powseries):

> sol:=powsolve({ODE,initvals});

sol := proc(powparm) ... end proc

y0x???O(x1516> tpsform(sol,x,16);

441641636912y0???y0x???y0x???y0x???y0x???98165615904913286025必须用_k, 这是powsolve使用的临时变量. > sol(_k);

?4a(_k???3)_k2)

也可以用powsolve给出的函数直接获得用递归形式定义的幂级数系数, 不过参数

2

例：求解一维谐振子的解：y???(??x)y?0 > alias(y=y(x)):

> ODE:=diff(y,x$2)+(epsilon-x^2)*y=0;

???????(????x2)y???0ODE := ?y2??x???2

> H:=powsolve(ODE);

H := proc(powparm) ... end proc

111??12?????30???24?C0???12C0?????60?C0??????> tpsform(H,x,8);

C0???C1x???12?C0x2???11?12?4?C1x3?????C0???C0???x???62412??1?1?52??120?C1???20C1??x?????111??1??728x6??????42???120?C1???20C1?????252?C1??x???O(x)????

> H(_k);

- 111 -

??a(_k???2)???a(_k???4)_k(_k???1)

2.5 常微分方程的数值解法

在对微分方程的解析解失效后, 可以求助于数值方法求解微分方程. 数值求解的好处是只要微分方程的条件足够多时一般都可求得结果, 然而所得结果是否正确则必须依赖相关数学基础加以判断. 调用函数dsolve求常微分方程初值问题的数值解时需加入参数type=numeric.

另一方面, 常微分方程初值问题数值求解还可以选择算法, 加入参数“method=方法参数”即可, 方法参数主要有：

rkf45：4~5阶变步长Runge-Kutta-Fehlberg法

dverk78：7~8阶变步长Runge-Kutta-Fehlberg法

classical：经典方法, 包括向前欧拉法, 改进欧拉法, 2、3、4阶龙格库塔法,

Sdams-Bashford方法等 gear：吉尔单步法 mgear：吉尔多步法 2.5.1变步长龙格库塔法

下面用4~5阶Runge-Kutta-Fehlberg法求解van der Pol方程：

?y???(1?y2)y??y?0 ??y(0)?0,y(0)??0.1?> ODE:=diff(y(t),t$2)-(1-y(t)^2)*diff(y(t),t)+y(t)=0;

????????y(t)???0????(1???y(t)2)?ODE := ?y(t)???ty(t)????t2?????2

> initvals:=y(0)=0,D(y)(0)=-0.1;

initvals := y(0)???0,D(y)(0)???-.1 > F:=dsolve({ODE,initvals},y(t),type=numeric);

F := proc(rkf45_x) ... end proc> F(0);

?t???0.,y(t)???0.,?y(t)???-.1???t???? ?

此时, 函数返回的是一个函数, 可以在给定的数值点上对它求值：

> F(1);

??t???1.,y(t)???-.144768589749425608?,y(t)???-.178104066128215944???? ?t?? 可以看到, F给出的是一个包括t、y(t)、D(y)(t)在内的有序表, 它对于每一个时间点

可以给出一组数值表达式. 有序表的每一项是一个等式, 可对其作图描述. > plot('rhs(F(t)[2])', t=0..15, title=\

- 112 -

> plots[odeplot](F,[t,y(t)],0..15,title=\

2.5.2吉尔法求解刚性方程

在科学和工程计算中, 常常会遇到这样一类常微分方程问题, 它可以表示成方程组：y??f(t,y),y(t0)?y0, 称其为刚性方程, 其解的分量数量相差很大, 分量的变化速度也相差很大. 如果用常规方法求解, 为了使变量有足够高的精度, 必须取很小的步长, 而为了使慢变分量达到近似的稳态解, 则需要很长的时间, 这样用小步长大时间跨度的计算, 必定造成庞大的计算量, 而且会使误差不断积累. 吉尔法是专门用来求解刚性方程的一种数值方法.

> ODE:=diff(u(t),t)=-2000*u(t)+999.75*v(t)+1000.25,diff(v(t),t)=u(t)-v(t);

??ODE := u(t)????2000u(t)???999.75v(t)???1000.25,v(t)???u(t)???v(t)

?t?t> initvals:=u(0)=0,v(0)=-2;

initvals := u(0)???0,v(0)???-2

> ansl:=dsolve({ODE,initvals},{u(t),v(t)},type=numeric,method=gear);

ansl := proc(x_gear) ... end proc

]

> ansl(10,0);

[t???10.,u(t)???.989893921726687442,v(t)???.979787842765888594> p1:=plots[odeplot] (ansl,[t,u(t)],0..20,color=red):

p2:=plots[odeplot] (ansl,[t,v(t)],0..20,color=blue):

plots[display] ({p1,p2}, title=\

- 113 -

2.5.3 经典数值方法

Maple中常微分方程数值解法中有一类被称作是“经典”(classical)方法. 当然, 称其为经典方法不是因为它们常用或是精度高, 而是因为它们的形式简单, 经常被用于计算方法课上的教学内容. 它们是一些常见的固定步长方法, 在dsolve中用参数method=classical[方法名称], 如果不特别指出, 将默认采用向前欧拉法. 主要有：

foreuler：向前欧拉法(默认)

hunform：Heun公式法(梯形方法, 改进欧拉法) imply：改进多项式法 rk2：二阶龙格库塔法 rk3：三阶龙格库塔法 rk4：四阶龙格库塔法

adambash：Adams-Bashford方法(预测法)

abmoulton：Adams-Bashford-Moulton方法(预测法) 下面给出微分方程数值方法的参数表：

参数名 参数类型 浮点数的一维数组 正整数 参数用途 指定初值向量 指定向量个数 Procedurelis：单个函数, 返回有序表 Listprocedure：函数的有序表 参数1：未知函数的个数 参数2：自变量 参数3：函数向量 参数4：导函数向量 参数用法 initial number output 'procedurelist'(默认) 指定生成单个函数或'listprocedure' 或多个函数的有序表 procedure 子程序名 用子程序形式指定第一尖常微分方程组的右边部分 start startinit value 浮点数 自变量起始值 布尔量(默认FALSE) 指定数值积分是否对dverk78不适用 总是从起始值开始 浮点数向量(一维数指定需要输出函数如果给定, 结果是一个2?2的矩阵. 元素[1,1]组) 值的自变量数值点 是一个向量, 含自变量名和函数名称; 元素[2,1]是一个数值矩阵, 其中第一列value的输入相同, 其他列中是相应的函数值 - 114 -

另外, 还有一些特殊的附加参数：

maxfun：整数类型, 用于最大的函数值数量, 默认值50000, 为负数时表示无限制 corrections：正整数类型, 指定每步修正值数量, 在abmoulton中使用, 建议值≤4 stepsize：浮点数值, 指定步长 下面看一个简单的例子：

> ODE:=diff(y(x),x)=y(x)-2*x/y(x);

?2xODE := y(x)???y(x)???

?xy(x)> initvals:=y(0)=1;

initvals := y(0)???1

> sol1:=dsolve({ODE,initvals},y(x),numeric,method=classical,stepsize=0.1,start=0);

sol1 := proc(x_classical) ... end proc

而其解析解为：

> sol2:=dsolve({diff(y(x),x)=y(x)-2*x/y(x), y(0)=1}, y(x));

sol2 := y(x)???2x???1

将两者图形同时绘制在同一坐标系中比较, 可以发现, 在经过一段时间后, 欧拉法的数值结果会产生较大误差.

> plot({rhs(sol2),'rhs(sol1(x)[2])'},x=0..2);

求解微方程, 无论使用什么方法或者加入什么选项, 求解完成后必须利用相关数学知识进行逻辑判断, 绝对不对简单迷信Maple给出的结果, 否则很有可能得到一个对于方程本身也许还看得过去, 但在数学或者物理意义上不合理的解.

2.6摄动法求解常微分方程

由于微分方程求解的复杂性, 一般微分方程常常不能求得精确解析解, 需要借助其它方法求得近似解或数值解, 或者两种方法兼而有之. 摄动法是重要的近似求解方法.

摄动法又称小参数法, 它处理含小参数?的系统, 一般当?=0时可求得解x0. 于是可把原系统的解展成?的幂级数x?x0?x1??x2?- 115 -

2??, 若这个级数当??0时一

致收敛，则称正则摄动, 否则称奇异摄动. 摄动法的种类繁多, 最有代表性的是庞加莱—林斯泰特(Poicare-Lindstedt)法, 在此, 我们以该方法求解van der Pol方程：

y????(1?y)y??y?0

2当?=0时该方程退化为数学单摆的常微分方程, 当?=1时为3.5讨论的情况, 对任意?, 该微分方程拥有一个渐进稳定的周期解, 称为极限环.

由于van der Pol方程中没有显式的时间依赖项, 不失一般性, 设初值为y(0)=0. 在庞加莱—林斯泰特法中, 时间通过变换拉伸：

????t, 其中?????ii?0i

对于y(?), van der Pol方程变为： ?y?????(1?y)y??y?0

22restart:

diff(y(t),t$2)-epsilon*(1-y(t)^2)*diff(y(t),t)+y(t)=0;

?????????(1???y(t)2)y(t)??t2???2??y(t)????y(t)???0 ????t???> ODE:=DEtools[Dchangevar]({t=tau/omega,y(t)=y(tau)},%,t,tau);

ODE := ?2????????y(?)???0??????(1???y(?)2)??y(?) ????y(?)?????2?????2> e_order:=6:

> macro(e=epsilon,t=tau):

> alias(seq(y[i]=eta[i](tau),i=0..e_order)): > e:=()->e:

> for i from 0 to e_order do

eta[i]:=t->eta[i](t) od:

> omega:=1+sum('w[i]*e^i','i'=1..e_order);

? := 1???w1????w2????w3????w4????w5????w6?

23456> y:=sum('eta[i]*e^i','i'=0..e_order);

y := ?0????1?????2?????3?????4?????5?????6?

23456> deqn:=simplify(collect(ODE,e),{e^(e_order+1)=0}):

- 116 -

> for i from 0 to e_order do ode[i]:=coeff(lhs(deqn),e,i)=0 od:

> ode[0];

???y0????2y0????0 ??????2> ode[1];

?????????????y???2w??y2???0 ????????????yyyy?20?11????0??0???21????0???????????22> ode[2];

2??????????y0??w1y0???2????y0??y0y1??????????????y1???????????y????y2??????22???22???????y1??y0??????????????y1????2?y0?w2???????y0??w1???2w1?22??????????????222????y?w???0???20?1??2

> dsolve({ode[0],eta[0](0)=0,D(eta[0])(0)=C[1]},eta[0](t));

y0???C1sin(?)

> eta[0]:=unapply(rhs(%),t);

?0 := ????C1sin(?)

> ode[1];

3????????Ccos(?)???y???2wCsin(?)???Ccos(?)sin(?)2???0 y11111???21???2> map(combine,ode[1],'trig');

3311????????Ccos(?)???y???2wCsin(?)???Ccos(?)???Ccos(3?)???0 y1111???21?4141??2> ode[1]:=map(collect,%,[sin(t),cos(t)]);

3?311?????ode1 := ?2w1C1sin(?)?????C???Ccos(?)???y???y???Ccos(3?)???0 ???11?11???21?44????2> dsolve({ode[1],eta[1](0)=0,D(eta[1])(0)=C[2]},eta[1](t),method=laplace);

1?sin(?)???y1???????8C1(C1???2)(C1???2)????w1C1???C2?????1C3???C?w?cos(?)???1C3cos(3?)???32111?321??

> map(collect,%,[sin(t),cos(t),t]);

1?sin(?)???y1???????8C1(C1???2)(C1???2)????w1C1???C2?????1C3???C?w?cos(?)???1C3cos(3?)???32111?321??

- 117 -

> solve({coeff(lhs(ode[1]),sin(t))=0,coeff(lhs(ode[1]),cos(t))=0});

{C1???0,w1???w1},{C1???2,w1???0},{C1???-2,w1???0}

> w[1]:=0:C[1]:=-2: > ode[1];

????????y???2cos(3?)???0 y1???21???2> dsolve({ode[1],eta[1](0)=0,D(eta[1])(0)=C[2]},eta[1](t),method=laplace); 11y1???cos(3?)???cos(?)???C2sin(?)

44> eta[1]:=unapply(rhs(%),tau); 11?1 := ????cos(3?)???cos(?)???C2sin(?)

44> map(combine,ode[2],'trig'):

> ode[2]:=map(collect,%,[sin(t),sin(3*t),cos(t),cos(3*t)]);

1?sin(?)???5sin(5?)???ode2 := ???4???4w2??4??2???3?2y2????sin(3?)???2C2cos(?)???3C2cos(3?)???y2???0????2??

> solve({coeff(lhs(ode[2]),sin(t))=0,coeff(lhs(ode[2]),cos(t))=0});

-1{C2???0,w2???}

16> assign(%):

> dsolve({ode[2],eta[2](0)=0,D(eta[2])(0)=C[3]},eta[2](t),method=laplace): > eta[2]:=unapply(rhs(%),t);

329?sin(?)???5sin(5?)?2 := ?????sin(3?)???????C ??963??1696??> for i from 0 to e_order do map(combine,ode[i],'trig'):

ode[i]:=map(collect,%,[seq(sin((2*j+1)*t),j=0..i),seq(cos((2*j+1)*t),j=0..i)]): solve({coeff(lhs(ode[i]),sin(t))=0,coeff(lhs(ode[i]),cos(t))=0}):

assign(%):

dsolve({ode[i],eta[i](0)=0,D(eta[i])(0)=C[i+1]},eta[i](t),method=laplace): collect(%,[seq(sin((2*j+1)*t),j=0..i),seq(cos((2*j+1)*t),j=0..i)]): eta[i]:=unapply(rhs(%),t); od: > omega;

1???116????2173072????435884736?

6- 118 -

> y(t):

> y:=unapply(simplify(y(t),{e^e_order=0}),t): > e:=1:y(t);

?1037927552960???sin(?)???12123732592001519540sin(?)cos(?)6???114941288008132sin(?)cos(?)4???6180sin(?)cos(?)8???13256sin(?)cos(?)2???52579573317760cos(?)3???55337200cos(?)11???912187129600cos(?)7???7999911105920cos(?)cos(?)5???cos(?)9

> plot(y(t),t=0..15);

在该问题的求解过程中, 前半部分我们按照交互式命令方式输入, 也就是把数学逻辑推理的过程“翻译”成Maple函数, 而在后半部分, 则采用程序设计方式书写了数学推导过程, 这是应用Maple解决实际问题的两种方式. 前一种方法只需了解Maple函数即可应用, 而后一种程序设计方式则需掌握Maple程序设计语言. 但是, 不论是那一种方式, 数学基础总是最重要的.

3 偏微分方程求解初步

Maple中偏微分方程求解器为pdsolve, 该函数及其它偏微分方程求解工具存于软件包PDEtools中. 函数pdsolve能够很快的辨认出偏微分方程是否为可以用标准方法求解的类型, 如果无法判别, 则pdsolve采用一种启发式的算法尝试偏微分方程按特征结构分离出来. pdsolve的策略就是寻找给定偏微分方程的通解, 如不能成功则寻找可以完全分离的变量, 因此, 该函数返回的结果可能为：

(i) 通解；

(ii) 近似的通解(即包含任意函数但又不足以得到通解的解)； (iii) 变量分离的非耦合的常微分方程.

如不能完全分离变量, 则函数再次调用自身, 如仍失败则返回未完全分离的变量并给出一个警告信息. 其命令格式为： pdsolve(PDE, f);

其中, PDE为偏微分方程, f为被求解函数.

下面通过几个例子说明pdsolve的用法： > PDE1:=diff(u(x,t),[t$2])-1/v^2*diff(u(x,t),[x$2]);

- 119 -

?222?x???PDE1 := ?2u(x,t)??????t???u(x,t)v2

> pdsolve(PDE1,u(x,t));

u(x,t)???_F1(?t???vx)???_F2(t???vx) > restart:

> PDE2:=a*x*diff(u(x,t),x)+b*t*diff(u(x,t),t)=0;

?????bt??u(x,t)????0PDE2 := ax?u(x,t) ??????x???t?????> pdsolve(PDE2,u(x,t));

t? u(x,t)???_F1??b??????a??????x?????> PDE3:=diff(u(x,t),x$2)-t^2*diff(u(x,t),t$2)-t*diff(u(x,t),t)=0;

???????t2PDE3 := ?u(x,t)2??x???2????????tu(x,t)??t2???2??u(x,t)????0 ????t???> pdsolve(PDE3,u(x,t));

u(x,t)???_F1(te)???_F2(tex(?x))

> restart:

> heatPDE:=diff(u(x,t),t)=diff(u(x,t),[x$2]);

heatPDE := ??tu(x,t)????22?xu(x,t)

> pdsolve(heatPDE,u(x,t));

2???(u(x,t)???_F1(x)_F2(t))&where?{_F1(x)???_c1_F1(x),_F2(t)???_c1_F2(t)}2??x?t?????

> dsolve(diff(F2(t),t)=c[1]*F2(t),F2(t));

(cF2(t)???_C1e1t)

(?cx)> dsolve(diff(F1(x),[x$2])=c[1]*F1(x));

(cF1(x)???_C1e1x)???_C2e1

(ct)> result:=rhs(%)*rhs(%%);

(cresult := (_C1e1x)(?c???_C2e- 120 -

1x))_C1e1

> subs(u(x,t)=result,heatPDE):expand(%);

(c(c21t)(c1x)_C1c1ee???_C1c1e(c1t)_C2???e(ct)(cx)1x)(c_C12c1e1e1???_C1c1e1t)_C2

> lhs(%)-rhs(%);

(cx)e10

- 121 -

本文来源：https://www.bwwdw.com/article/l8qx.html