也可以用powsolve给出的函数直接获得用递归形式定义的幂级数系数, 不过参数必须用_k, 这是powsolve使用的临时变量. > sol(_k); a(_k???3)?42_k 2??例：求解一维谐振子的解：y?(??x)y?0 > alias(y=y(x)): > ODE:=diff(y,x$2)+(epsilon-x^2)*y=0; 2???2?ODE := ?2y????(????x)y???0? ?x??> H:=powsolve(ODE); H := proc(powparm) ... end proc- 125 -
> tpsform(H,x,8); 111211211?1211???????2345C0???C1x????C0x????C1x?????C0???C0x????C1???C1x??????C0???C0????C0????????24??120??30?24?60?26122012????????1?121????1?C1?x7???O(x8) x6???????C1???C1????42??120??20???252?> H(_k); ?a(_k???2)???a(_k???4)? _k(_k???1)2.5 常微分方程的数值解法
在对微分方程的解析解失效后, 可以求助于数值方法求解微分方程. 数值求解的好处是只要微分方程的条件足够多时一般都可求得结果, 然而所得结果是否正确则必须依赖相关数学基础加以判断. 调用函数dsolve求常微分方程初值问题的数值解时需加入参数type=numeric.
另一方面, 常微分方程初值问题数值求解还可以选择算法, 加入参数“method=方法参数”即可, 方法参数主要有： rkf45：4~5阶变步长Runge-Kutta-Fehlberg法
- 126 -
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; 2???2ODE := ?2y(t)?????(1???y(t))??t????y(t)????y(t)???0????t???> initvals:=y(0)=0,D(y)(0)=-0.1;
> F:=dsolve({ODE,initvals},y(t),type=numeric);
F := proc(rkf45_x) ... end proc
此时, 函数返回的是一个函数, 可以在给定的数值点上对它求值：
initvals := y(0)???0,D(y)(0)???-.1- 127 -
> 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=\ > plots[odeplot](F,[t,y(t)],0..15,title=\of the Van de Pol's Equation\
- 128 -
2.5.2吉尔法求解刚性方程 在科学和工程计算中, 常常会遇到这样一类常微分方程问题, 它可以表示成方程组：y??f(t,y),y(t)?y, 称其为刚性方程, 其解的分量数量相差很大, 分量的变化速度也相差很大. 如果用常规方法求解, 为了使变量有足够高的精度, 必须取很小的步长, 而为了使慢变分量达到近似的稳态解, 则需要很长的时间, 这样用小步长大时间跨度的计算, 必定造成庞大的计算量, 而且会使误差不断积累. 吉尔法是专门用来求解刚性方程的一种数值方法. > 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?t00 > initvals:=u(0)=0,v(0)=-2; initvals := u(0)???0,v(0)???-2 - 129 -