Maple笔记2--常微分方程求解

转需看

原文地址：Maple笔记2--常微分方程求解作者：Lionel

来源：网络论坛转载（VB资料库）

常微分方程求解
微分方程求解是数学研究与应用的一个重点和难点. Maple能够显式或隐式地解析地求解许多微分方程求解.在常微分方程求解器dsolve中使用了一些传统的技术例如laplace变换和积分因子法等,函数pdesolve则使用诸如特征根法等经典方法求解偏微分方程. 此外, Maple还提供了可作摄动解的所有工具,例如Poincare-Lindstedt法和高阶多重尺度法. 
帮助处理常微分方程(组)的各类函数存于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—依赖于要求解的问题类型. 
例如, 对于一阶常微分方程 可用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之前应将表达式化简, 否则可能会出错.
再看下面常微分方程的求解： 
> 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)} );

2.4 常微分方程的级数解法
1) 泰勒级数解法
当一个常微分方程的解析解难以求得时, 可以用Maple求得方程解的级数近似, 这在大多数情况下是一种非常好的方法.级数解法是一种半解析半数值的方法. 
泰勒级数法的使用命令为：
dsolve({ODE,Ics}, y(x), 'series'); 或dsolve({ODE,Ics}, y(x),'type=series');
下面求解物理摆的大幅振动方程： , 其中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中.但是, 这一求解函数的使用范围很有限,它只可以用来求解多项式系数的线性常微分方程或方程组,其求解命令为：powseries[function](prep)或直接载入软件包后用function(prep),prep为求解的线性微分方程及其初值. 
例：求解： 
>ODE:=x*diff(y(x),x$2)+diff(y(x),x)+4*x^2*y(x)=0;

> dsolve(ODE,y(x));

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

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

> tpsform(sol,x,16);

也可以用powsolve给出的函数直接获得用递归形式定义的幂级数系数, 不过参数必须用_k,这是powsolve使用的临时变量. 
> sol(_k);

例：求解一维谐振子的解： 
> alias(y=y(x)):
> ODE:=diff(y,x$2)+(epsilon-x^2)*y=0;

> H:=powsolve(ODE);

> tpsform(H,x,8);

> H(_k);

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方程：

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

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

> F:=dsolve({ODE,initvals},y(t),type=numeric);

    此时, 函数返回的是一个函数,可以在给定的数值点上对它求值：
> F(0);

> F(1);

    可以看到,F给出的是一个包括t、y(t)、D(y)(t)在内的有序表, 它对于每一个时间点可以给出一组数值表达式. 有序表的每一项是一个等式,可对其作图描述. 
> plot('rhs(F(t)[2])', t=0..15, title="solution of the Van dePol's Equation");

> plots[odeplot](F,[t,y(t)],0..15,title="solution of the Van dePol's Equation");

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

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

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

> ansl(10,0);

> 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="Solution of a stiffequation");

2.5.3 经典数值方法
Maple中常微分方程数值解法中有一类被称作是“经典”(classical)方法. 当然, 称其为经典方法不是因为它们常用或是精度高,而是因为它们的形式简单, 经常被用于计算方法课上的教学内容. 它们是一些常见的固定步长方法,在dsolve中用参数method=classical[方法名称], 如果不特别指出, 将默认采用向前欧拉法.主要有：
foreuler：向前欧拉法(默认)
hunform：Heun公式法(梯形方法, 改进欧拉法)
imply：改进多项式法
rk2：二阶龙格库塔法
rk3：三阶龙格库塔法
rk4：四阶龙格库塔法
adambash：Adams-Bashford方法(预测法)
abmoulton：Adams-Bashford-Moulton方法(预测法)
下面给出微分方程数值方法的参数表：

参数名

参数类型

参数用途

参数用法

initial
浮点数的一维数组
指定初值向量

number
正整数
指定向量个数

output
'procedurelist'(默认)
或'listprocedure'
指定生成单个函数或多个函数的有序表
Procedurelis：单个函数, 返回有序表
Listprocedure：函数的有序表
procedure
子程序名
用子程序形式指定第一尖常微分方程组的右边部分
参数1：未知函数的个数
参数2：自变量
参数3：函数向量
参数4：导函数向量
start
浮点数
自变量起始值

startinit
布尔量(默认FALSE)
指定数值积分是否总是从起始值开始
对dverk78不适用
value
浮点数向量(一维数组)
指定需要输出函数值的自变量数值点
如果给定, 结果是一个 的矩阵. 元素[1,1]是一个向量, 含自变量名和函数名称; 元素[2,1]是一个数值矩阵, 其中第一列value的输入相同, 其他列中是相应的函数值


另外, 还有一些特殊的附加参数：
maxfun：整数类型, 用于最大的函数值数量, 默认值50000, 为负数时表示无限制
corrections：正整数类型, 指定每步修正值数量, 在abmoulton中使用, 建议值≤4
stepsize：浮点数值, 指定步长
下面看一个简单的例子：
> ODE:=diff(y(x),x)=y(x)-2*x/y(x);

> initvals:=y(0)=1;

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

    而其解析解为：
> sol2:=dsolve({diff(y(x),x)=y(x)-2*x/y(x), y(0)=1},y(x));

    将两者图形同时绘制在同一坐标系中比较, 可以发现,在经过一段时间后, 欧拉法的数值结果会产生较大误差. 
> plot({rhs(sol2),'rhs(sol1(x)[2])'},x=0..2);

    求解微方程, 无论使用什么方法或者加入什么选项,求解完成后必须利用相关数学知识进行逻辑判断, 绝对不对简单迷信Maple给出的结果, 否则很有可能得到一个对于方程本身也许还看得过去,但在数学或者物理意义上不合理的解. 
2.6摄动法求解常微分方程
由于微分方程求解的复杂性, 一般微分方程常常不能求得精确解析解, 需要借助其它方法求得近似解或数值解, 或者两种方法兼而有之.摄动法是重要的近似求解方法.  
摄动法又称小参数法, 它处理含小参数 的系统, 一般当 =0时可求得解x0. 于是可把原系统的解展成 的幂级数 ,若这个级数当  0时一致收敛，则称正则摄动, 否则称奇异摄动.摄动法的种类繁多, 最有代表性的是庞加莱—林斯泰特(Poicare-Lindstedt)法, 在此, 我们以该方法求解van derPol方程：

当 =0时该方程退化为数学单摆的常微分方程, 当 =1时为3.5讨论的情况, 对任意 , 该微分方程拥有一个渐进稳定的周期解,称为极限环. 
由于van der Pol方程中没有显式的时间依赖项, 不失一般性, 设初值为y(0)=0. 在庞加莱—林斯泰特法中,时间通过变换拉伸：
, 其中 
对于 , van der Pol方程变为：

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

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

> e_order:=6:
> macro(e=epsilon,t=tau):
> alias(seq(y=eta(tau),i=0..e_order)):
> e:=()->e:
> for i from 0 to e_order do
eta:=t->eta(t)
od:
> omega:=1+sum('w*e^i','i'=1..e_order);

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

> deqn:=simplify(collect(ODE,e),{e^(e_order+1)=0}):
> for i from 0 to e_order do
ode:=coeff(lhs(deqn),e,i)=0
od:
> ode[0];

> ode[1];

> ode[2];

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

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

> ode[1];

> map(combine,ode[1],'trig');

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

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

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

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

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

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

> eta[1]:=unapply(rhs(%),tau);

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

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

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

> for i from 0 to e_order do
  map(combine,ode,'trig'):
  ode:=map(collect,%,[seq(sin((2*j+1)*t),j=0..i),seq(cos((2*j+1)*t),j=0..i)]):
  solve({coeff(lhs(ode),sin(t))=0,coeff(lhs(ode),cos(t))=0}):
   assign(%):
  dsolve({ode,eta(0)=0,D(eta)(0)=C[i+1]},eta(t),method=laplace):
  collect(%,[seq(sin((2*j+1)*t),j=0..i),seq(cos((2*j+1)*t),j=0..i)]):
  eta:=unapply(rhs(%),t);
od:
> omega;

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

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

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