第四章 方程求解(2)

帮助处理常微分方程(组)的各类函数存于Detools软件包中，函数种类主要有：可视化类的函数，处理宠加莱动态系统的函数，调整微分方程的函数，处理积分因子、李对称法和常微分方程分类的函数，微分算子的函数，利用可积性与微分消去的方法简化微分方程的函数，以及构造封闭解的函数等。更重要的是其提供的强大的图形绘制命令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'?ylnxy?y(可用dsolve直接求得解析解： > ODE:=x*diff(y(x),x)=y(x)*ln(x*y(x))-y(x);

> dsolve(ODE,y(x));

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

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

> dsolve(ODE,y);

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

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

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

> evalb(%);

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

再看下面常微分方程的求解：12+=′yy > alias(y=y(x)):

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

> dsolve(ODE,y);

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

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

> sol:=dsolve(ODE,y);

由此可见，对于不能表示成显式结果的微分方程解，Maple尽可能将结果表示成隐式解。另外，对于平凡解y=0常常忽略，这一点应该引起注意。

dsolve对于求解微分方程初值问题也十分方便的: > ODE:=diff(u(t),t$2)+omega^2*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)+omega^2*u(t)=Heaviside(t);

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

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

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

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

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

> with(plots):

plot3d(rhs(%%),omega=2/3..4/3,t=0..20,style=hidden,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}, {vars}); 其中, 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)=t^2+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()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)=v[0]/l;

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

> Order:=11:

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

2) 幂级数解法

对于一个符号代数系统来说，幂级数是必不可少的微分方程求解工具。幂级数求解函数powsolve存于工具包powseries中，但是，这一求解函数的使用范围很