摄动法又称小参数法, 它处理含小参数?的系统, 一般当?=0时可求得解x0. 于是可把原系统的解展成?的幂级数x?x?x??x???, 若这个级数当??0时一致收敛，则称正则摄动, 否则称奇异摄动. 摄动法的种类繁多, 最有代表性的是庞加莱—林斯泰特(Poicare-Lindstedt)法, 在此, 我们以该方法求解van der Pol方程：
y????(1?y)y??y?0
20122当?=0时该方程退化为数学单摆的常微分方程, 当?=1时为3.5讨论的情况, 对任意?, 该微分方程拥有一个渐进稳定的周期解, 称为极限环.
由于van der Pol方程中没有显式的时间依赖项, 不失一般性, 设初值为y(0)=0. 在庞加莱—林斯泰特法中, 时间通过变换拉伸：
???t, 其中?????
?iii?0对于y(?), van der Pol方程变为： ?y?????(1?y)y??y?0
22restart:
- 135 -
diff(y(t),t$2)-epsilon*(1-y(t)^2)*diff(y(t),t)+y(t)=0; 2????????y(t)???0?2y(t)?????(1???y(t)2)?y(t)????t???t????? > ODE:=DEtools[Dchangevar]({t=tau/omega,y(t)=y(tau)},%,t,tau); 2????????(1???y(?)2)?ODE := ??y(?)???2???2??y(?)????y(?)???0????????> e_order:=6:
> macro(e=epsilon,t=tau):
> alias(seq(y[i]=eta[i](tau),i=0..e_order)): > e:=()->e:
> for i from 0 to e_order do eta[i]:=t->eta[i](t) od:
> omega:=1+sum('w[i]*e^i','i'=1..e_order);
? := 1???w1????w2?2???w3?3???w4?4???w5?5???w6?6
> y:=sum('eta[i]*e^i','i'=0..e_order);
- 136 -
y := ?0????1?????2?2????3?3????4?4????5?5????6?6
> deqn:=simplify(collect(ODE,e),{e^(e_order+1)=0}): > for i from 0 to e_order do ode[i]:=coeff(lhs(deqn),e,i)=0 od: > ode[0]; 2????y0????y???20????0?? > ode[1];
22????????????2?2y1???????????y???y???2wy?20?0?11????y0??y0???0?????????????????> ode[2];
?2??2??????y?wy2???2??y?yy?????y?????????y2??????????????0?10???0?01???1????2y2????y1??y0???????y0??w1???2w1?????????????222?????????2?2y1????2?2y0?w2????2y0?w1???0??????????????????> dsolve({ode[0],eta[0](0)=0,D(eta[0])(0)=C[1]},eta[0](t));
y???Csin(?)
01> eta[0]:=unapply(rhs(%),t);
- 137 -
?0 := ????C1sin(?)
> ode[1];
23????2y1????C1cos(?)???y1???2w1C1sin(?)???C1cos(?)sin(?)2???0??????> map(combine,ode[1],'trig'); 21313????2y1????C1cos(?)???y1???2w1C1sin(?)???C1cos(?)???C1cos(3?)???0????44??> ode[1]:=map(collect,%,[sin(t),cos(t)]); 213?13??????ode1 := ?2w1C1sin(?)?????C???Ccos(?)???y???y???C1cos(3?)???021??141??1?4??????> dsolve({ode[1],eta[1](0)=0,D(eta[1])(0)=C[2]},eta[1](t),method=laplace); 1?sin(?)????1C3???C?w?cos(?)???1C3cos(3?)y1?????C(C???2)(C???2)????wC???C???8111112???32111??321????> map(collect,%,[sin(t),cos(t),t]); 1?sin(?)????1C3???C?w?cos(?)???1C3cos(3?)y1????????8C1(C1???2)(C1???2)????w1C1???C2???32111??321???? 111> solve({coeff(lhs(ode[1]),sin(t))=0,coeff(lhs(ode[1]),cos(t))=0});
{C???0,w???w},{C???2,w???0},{C???-2,w???0}
1111- 138 -
> w[1]:=0:C[1]:=-2: > ode[1];
2????2y1????y1???2cos(3?)???0?????? > dsolve({ode[1],eta[1](0)=0,D(eta[1])(0)=C[2]},eta[1](t),method=laplace); 11y???cos(3?)???cos(?)???Csin(?) 4412> eta[1]:=unapply(rhs(%),tau); 11?1 := ????cos(3?)???cos(?)???C2sin(?)44 > map(combine,ode[2],'trig'): > ode[2]:=map(collect,%,[sin(t),sin(3*t),cos(t),cos(3*t)]); 1?2?3?sin(?)???5sin(5?)??????ode2 := ????4w??42?????2y2????2sin(3?)???2C2cos(?)???3C2cos(3?)???y2???04????> solve({coeff(lhs(ode[2]),sin(t))=0,coeff(lhs(ode[2]),cos(t))=0}); -1{C???0,w???} 1622> assign(%): > dsolve({ode[2],eta[2](0)=0,D(eta[2])(0)=C[3]},eta[2](t),method=laplace): - 139 -