B3-Maple中的偏微分方程求解

B3: Maple中的偏微分方程求解 西希安工程模拟软件（上海）有限公司，2008 11.0 Maple中的微分方程求解器介绍 Maple中微分方程求解器使用领先的算法求解以下问题： 常微分方程 (ODEs): dsolve 命令用于求解线性和非线性ODEs, 初始值问题 (IVP), 以及边界 值问题 (BVP)，可以通过参数项选择求符号解 (解析解) 或数值解。ODE Analyzer Assistant 微分方程分析器助手提供一个交互式用户界面方便用户求解 ODE 以及显示结果的图形。了 解更多信息，参考帮助系统中的 dsolve, dsolve/numeric, 和 ODE Analyzer. 偏微分方程 (PDEs): pdsolve 命令用于求 PDEs 和含边界值问题的 PDEs 的符号解或数值 解。使用Maple的PDE 工具可以完成对PDE 系统的结构分析和指数降阶处理。了解更多信 息，参考帮助系统中的 pdsolve and pdsolve/numeric. 微分-代数方程 (DAEs): dsolve/numeric 命令是符号-数值混合求解器，使用符号预处理和降 阶技术，让Maple能够求解高指数的DAE问题。Maple内置三个求解器用于处理DAEs：1） 修正的 Runge-Kutta Fehlberg 方法，2）Rosenbrock 方法，以及 3）修正的拓展后向差分隐 式方法。 11.1 求解偏微分方程PDE 问题(BVP和IVP) Maple 求解经典力学难题的能力是非常著名的，它的数值和符号偏微分方程求解器是其中的 重要工具。 例子：在不同的边界条件下，求波动方程的数值解、解析解、和图形解。 11.1.1 初始化 下面的Maple代码定义了一个名为PX的程序，生成函数的周期展开。 PX := proc(h::{algebraic,procedure},g::{range,name=range}) local L, D, var; if type(g,'range') then L := lhs(g); D := rhs(g) - L; (1)(1) (2)(2) if not type(h,'procedure') then var := indets(h,'name'); if nops(var) <> 1 then error "need to specify a variable"; end if; var := op(var); end if; else L := lhs(rhs(g));D := rhs(rhs(g)) - L;var := lhs(g); end if; if type(h,'procedure') then proc(x::algebraic) h(x - floor((x-L)/D)*D); end; else proc(x::algebraic) eval(h, var = x - floor((x-L)/D)*D); end; end if; end: 11.1.2 数值解和图形解 一个空间变量的波动方程是： 假设初始形状由下面的函数给出： 对应的图形如下： (4)(4) (3)(3) x 0 2 4 6 8 10 0 假设波动方程的初始条件和边界条件如下： Maple中的命令pdsolve将求解单变量演化方程（双曲和抛物）的数值解（有限差分）。 这个命令创建了一个模块，可以看到模块的输出函数是 plot, plot3d, animate, 和 value。 时间 的图形： x 0 2 4 6 8 10 0 图 2 解的表面图： Figure 3 动画： (5)(5) (6)(6) x 2 4 6 8 10 0 1 time = 0.000000 给出计算值 的函数： 例如， 的解： 0.364667328524451 11.1.3 D'Alembert 解析解 对于这个问题，D'Alembert 解的形式是： = (7)(7) (8)(8) 因此 的精确值是： 0.3750000000 在 区间 上的奇次展开如下： 通过【图4】验证。 (9)(9) x 5 10 图 4 的奇次周期展开如下： 我们注意到输出是一个程序体，忽略Maple输出的细节，我们再次使用图形来验证工作， 得到【图5】。 x 0 10 20 30 图 5 最后，通过下式组合解 的输出是一个分段函数的复杂组合，但是下面的动画显示了在指定区间内波的相应运 动。 x 5 10 u 1 【图6】是解的表面，与动画提供的信息一致。 (11)(11) (10)(10) 图 6 11.1.4 扩展（可选）：传播和Klein-Gordon方程 10.2.4.1 数值解 波和Klein-Gordon方程分别表示为： 两者之间的区别是Klein-Gordon方程中有 项。 (12)(12) (13)(13) 为了说明两个方程解之间的查表，考虑下面的初始扰动： 传播长度设定为2 0，端点固定。初始条件和边界条件的方程是： 方程的数值解分别是： 方程的动画分别是： 波动方程的动画是： x 10 20 30 40 0 1 Klein-Gordon传播方程的动画是： (15)(15) (14)(14) x 10 20 30 40 0 1 10.2.4.2 解析解 Klein-Gordon方程产生分离的变量。因此，假设解的形式如下： 因此Klein-Gordon方程，这里 ，变为： 除以 得到： (20)(20) (21)(21) (15)(15) (17)(17) (16)(16) (19)(19) (18)(18) -1 移到左边： 基于以前的经验，我们引入Bernoulli分离常数 ，这两个常微分方程变为： 在不同的假设下，得到边界条件： = 0 使用其中的第一个得到： 使用其中的第二个得到： 第二个常微分方程变为： (21)(21) (15)(15) (26)(26) (27)(27) x 10 20 30 40 0 1 可以发现数值解和级数近似值几乎等同。 11.1.5 求解初始值偏微分方程问题（热传导方程） 数值解 pde := diff(u(x, t), t) = diff(u(x, t), x, x); 初始条件： ibc := [(D[1](u))(0, t)+a*u(0, t) = h(t), u(x, 0) = 0]; (21)(21) (28)(28) (29)(29) (15)(15) 可以得到数值解，但是我们需要定义 a 和 h 的值，以及提供第二个边界条件： ibc1:=[D[1](u)(0,t)+2*u(0,t)=cos(t),u(x,0)=0,u(0,t)=1]; pds:=pdsolve(pde,ibc1,numeric,time=t,range=0..1); 我们可以做其他计算，例如生成解的三维图形： pds:-plot3d(t=0..1,x=0..1,axes=boxed); 关于偏微分方程组数值解的更多信息，请参阅下面的帮助文件： ?pdsolve,numeric 上面的帮助页中包含了处理热传导方程的示例。 解析解 (21)(21) (33)(33) (15)(15) (34)(34) (30)(30) (32)(32) (31)(31) 我们可以使用 pdsolve （无需参数项）命令计算偏微分方程的解析解： gsol:=pdsolve(pde); 结果可以拆分为方程（eq）和条件（conds）： eq:=op(gsol)[1]; conds:=op(op(gsol)[2]); 我们可以用 dsolve 命令求解常微分方程组，然后使用结果对 eq 求值： conditions:=dsolve(conds); solution:=eval(eq,conditions); 我们可以判定结果是否满足偏微分方程： pdetest(solution,pde); 0 然后我们可以调整常数 _c1, _C1, _C2 和 _C3 满足期望的条件。 关于 pdsolve 命令的更多信息，请参考下面的帮助页： ?pdsolve 其他帮助 pdetest, pdsolve