第四章 方程求解(3)

有限，它只可以用来求解多项式系数的线性常微分方程或方程组，其求解命令为：

powseries[function] (prep)

或直接载入软件包后用function(prep), prep为求解的线性微分方程及其初值. 例：求解： 042=+′′+′yxyyx

> 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=\

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

2.5.2吉尔法求解刚性方程

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

> 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);

:= anslproc()... end procx_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=\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 指定生Procedurelis：单个函数, 返回'procedurelist'(默认) 成单个有序表 函数或Listprocedure：函数的有序表 或'listprocedure' 多个函数的有序表 procedure 子程序名 用子程参数1：未知函数的个数 序形式参数2：自变量 指定第参数3：函数向量 一尖常参数4：导函数向量 微分方程组的右边部分 start 浮点数 自变量起始值 startinit 布尔量(默认FALSE) 指定对dverk78不适用 数值积分是否总是从起始值开始 value 浮点数向量(一维数指定如果给定, 结果是一个22×的矩组) 需要阵. 元素[1,1]是一个向量, 含输出自变量名和函数名称; 元素函数[2,1]是一个数值矩阵, 其中第值的一列value的输入相同, 其他列自变中是相应的函数值 量数值点 另外, 还有一些特殊的附加参数： maxfun：整数类型, 用于最大的函数值数量, 默认值50000, 为负数时表示无限制

corrections：正整数类型, 指定每步修正值数量, 在abmoulton中使用, 建议值≤4

stepsize：浮点数值, 指定步长 下面看一个简单的例子： > ODE:=diff(y(x),x)=y(x)-2*x/y(x);

:= ODE = ??x()yx ? ()yx2x()yx

> initvals:=y(0)=1;

:= initvals = ()y01

>

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. 于是可把原系统的解展成ε的幂级数, 若这个级数当?????2210εεxxxxε→0时一致收敛，则称正则摄动, 否则称奇异摄动. 摄动法的种类繁多, 最有代表性的是庞加莱—林斯泰特(Poicare-Lindstedt)法, 在此, 我们以该方法求解van der Pol方程：

0)1(2=+′??′′yyyyε

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

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

tωτ=, 其中Σ∞==0iiiεωω 对于)(τy, van der Pol方程变为： 0)1(22=+′??′′yyyy …… 此处隐藏：2532字，全部文档内容请下载后查看。喜欢就下载吧 ……