1vandpol.mac is from the book "Computer Algebra in Applied Mathematics: 2An introduction to MACSYMA", by Richard H Rand, Pitman (1984). 3 4The van der Pol equation is 5 6 x''(T) + x - e*(1-x^2)*x'(T) = 0 7 8 9This program uses Lindstedt's method to find the limit cycle, the 10asymptotically stable periodic solution. Time is "stretched" using 11the transformation 12 t = w(e)*T 13where 14 w(e) = 1 + k[1]*e + k[2]*e^2 + ... 15We then expand x(t) as a power series in e 16 x(t) = y[0](t) + e*y[1] + e^2*y[2] + ... 17substitute this into the de and equate coefficients of e 18 19The original code produces an error. A small patch (below) is 20required to get the code to run with maxima-5.9.0cvs. (A similar patch 21is needed for duffing.mac). 22 23The run below, using maxima-5.9.0cvs, reproduces the result on pages 2481-82 of the book. 25 26(C1) load("./vandpol.mac"); 27(D1) ./vandpol.mac 28(C2) vanderpol(6); 29 3 SIN(t) SIN(3 t) 30y (t) = -------- - -------- 31 1 4 4 32 33 5 COS(5 t) 3 COS(3 t) COS(t) 34y (t) = - ---------- + ---------- - ------ 35 2 96 16 8 36 37 7 SIN(7 t) 35 SIN(5 t) 21 SIN(3 t) 7 SIN(t) 38y (t) = ---------- - ----------- + ----------- - -------- 39 3 576 576 256 256 40 41 61 COS(9 t) 2149 COS(7 t) 1085 COS(5 t) 47 COS(3 t) 73 COS(t) 42y (t) = ----------- - ------------- + ------------- - ----------- + --------- 43 4 20480 110592 27648 1536 12288 44 45 5533 SIN(11 t) 7457 SIN(9 t) 110621 SIN(7 t) 52885 SIN(5 t) 46y (t) = - -------------- + ------------- - --------------- + -------------- 47 5 7372800 1228800 6635520 2654208 48 49 2591 SIN(3 t) 12971 SIN(t) 50 - ------------- - ------------ 51 294912 4423680 52 53 6 4 2 54 35 e 17 e e 55w= ------ + ----- - -- + 1 56 884736 3072 16 57 58 59The required patch is 60 61--- vandpol.mac.orig 2004-02-14 11:48:48.478448000 +1100 62+++ vandpol.mac 2004-02-14 11:57:35.626449600 +1100 63@@ -43,7 +43,8 @@ 64 IF i = 1 THEN f[i]:solve(coeff(temp1,cos(t)),k[i]) 65 ELSE f[i]:solve([coeff(temp1,cos(t)),coeff(temp1,sin(t))], 66 [k[i],b[i-1]]),temp1:ev(temp1,f[i]))$ 67-step3(i):=temp1:ev(ode2(temp1,y[i](t),t),%k1:a[i],%k2:b[i])$ 68+step3(i):=(temp1:ode2(temp1,y[i](t),t), 69+ temp1:subst(a[i],%k1,temp1),temp1:subst(b[i],%k2,temp1))$ 70 step4(i):=(temp2:rhs(temp1),temp2:diff(temp2,t), 71 temp2:solve(ev(temp2,t:0),a[i]))$ 72 step5(i):=e[i]:ev(temp1,temp2)$ 73 74 75Local Variables: *** 76mode: Text *** 77End: ***