1duffing.mac is from the book "Computer Algebra in Applied Mathematics: 2An introduction to MACSYMA", by Richard H Rand, Pitman (1984). 3 4Duffing's equation is 5 6 x''(T) + x + e*x^3 = e*f*cos(w*T) 7 8It is a model of a nonlinear oscillator driven by a sinusiodal forcing 9function. When e is small and the forcing frequency is close to the 10natural frequency (of unity) the equation has one or more steady state 11periodic solutions. 12 13This program determines the periodic solutions using Lindstedt's method. 14The code is very similar to vandpol.mac and mathieu.mac. 15 16The original code produces an error. A small patch (below) is required 17to get the code to run with maxima-5.9.0cvs. 18 19The run below, using maxima-5.9.0cvs, reproduce the result on pages 20127-134 of the book. In the first part, the problem is solved step by 21step. Then the whole procedure is automated. 22 23(C1) load("./duffing.mac"); 24(D1) ./duffing.mac 25(C2) setup1(2); 26(D2) DONE 27(C3) w; 28 2 29(D3) k e + k e + 1 30 2 1 31(C4) x; 32 2 33(D4) a COS(t) + e y (t) + e y (t) 34 2 1 35(C5) setup2(2); 36(D5) DONE 37(C6) eq[1]; 38 2 39 d 3 3 40(D6)/R/ --- (y (t)) + a COS (t) + (- f - 2 k a) COS(t) + y (t) 41 2 1 1 1 42 dt 43(C7) eq[2]; 44 2 45 d 3 3 2 2 46(D7)/R/ --- (y (t)) - 2 k a COS (t) + 3 a y (t) COS (t) 47 2 2 1 1 48 dt 49 50 2 51 + (2 k f + (- 2 k + 3 k ) a) COS(t) + y (t) - 2 k y (t) 52 1 2 1 2 1 1 53(C8) step1(1); 54 2 3 3 55 d a COS(3 t) 3 a COS(t) 56(D8) --- (y (t)) + ----------- - f COS(t) + ----------- - 2 k a COS(t) + y (t) 57 2 1 4 4 1 1 58 dt 59(C9) step2(1); 60 2 3 3 3 61 d a COS(3 t) (4 f - 3 a ) COS(t) 3 a COS(t) 62(D9) --- (y (t)) + ----------- + ------------------- - f COS(t) + ----------- 63 2 1 4 4 4 64 dt 65 66 + y (t) 67 1 68(C10) f[1]; 69 3 70 4 f - 3 a 71(D10) [k = - ----------] 72 1 8 a 73(C11) step3(1); 74 3 75 a COS(3 t) 76(D11) y (t) = ----------- + %K1 SIN(t) + %K2 COS(t) 77 1 32 78(C12) step4(1); 79Inconsistent equations: (1 2) 80#0: step4(i=1) 81 -- an error. Quitting. To debug this try DEBUGMODE(TRUE);) 82(C13) 83 84 85After a small patch 86 87(C1) load("./duffing.mac"); 88(D1) ./duffing.mac 89(C2) setup1(2); 90(D2) DONE 91(C3) w; 92 2 93(D3) k e + k e + 1 94 2 1 95(C4) x; 96 2 97(D4) a COS(t) + e y (t) + e y (t) 98 2 1 99(C5) setup2(2); 100(D5) DONE 101(C6) eq[1]; 102 2 103 d 3 3 104(D6)/R/ --- (y (t)) + a COS (t) + (- f - 2 k a) COS(t) + y (t) 105 2 1 1 1 106 dt 107(C7) eq[2]; 108 2 109 d 3 3 2 2 110(D7)/R/ --- (y (t)) - 2 k a COS (t) + 3 a y (t) COS (t) 111 2 2 1 1 112 dt 113 114 2 115 + (2 k f + (- 2 k + 3 k ) a) COS(t) + y (t) - 2 k y (t) 116 1 2 1 2 1 1 117(C8) step1(1); 118 2 3 3 119 d a COS(3 t) 3 a COS(t) 120(D8) --- (y (t)) + ----------- - f COS(t) + ----------- - 2 k a COS(t) + y (t) 121 2 1 4 4 1 1 122 dt 123(C9) step2(1); 124 2 3 3 3 125 d a COS(3 t) (4 f - 3 a ) COS(t) 3 a COS(t) 126(D9) --- (y (t)) + ----------- + ------------------- - f COS(t) + ----------- 127 2 1 4 4 4 128 dt 129 130 + y (t) 131 1 132(C10) f[1]; 133 3 134 4 f - 3 a 135(D10) [k = - ----------] 136 1 8 a 137(C11) step3(1); 138 3 139 a COS(3 t) 140(D11) y (t) = ----------- + a SIN(t) + b COS(t) 141 1 32 1 1 142(C12) step4(1); 143 3 144 a 145(D12) [[a = 0, b = - --]] 146 1 1 32 147(C13) step5(1); 148 3 3 149 a COS(3 t) a COS(t) 150(D13) y (t) = ----------- - --------- 151 1 32 32 152(C14) step1(2); 153 2 5 2 5 2 154 d 3 a COS(5 t) 9 a f COS(3 t) 3 a COS(3 t) f COS(t) 155(D14) --- (y (t)) + ------------- + --------------- - ------------- - --------- 156 2 2 128 32 16 4 a 157 dt 158 159 2 5 160 11 a f COS(t) 21 a COS(t) 161 + -------------- - ------------ - 2 k a COS(t) + y (t) 162 32 128 2 2 163(C15) step2(2); 164 2 5 2 5 165 d 3 a COS(5 t) 9 a f COS(3 t) 3 a COS(3 t) 166(D15) --- (y (t)) + ------------- + --------------- - ------------- 167 2 2 128 32 16 168 dt 169 170 2 3 6 2 2 5 171 (32 f - 44 a f + 21 a ) COS(t) f COS(t) 11 a f COS(t) 21 a COS(t) 172 + -------------------------------- - --------- + -------------- - ------------ 173 128 a 4 a 32 128 174 175 + y (t) 176 2 177(C16) f[2]; 178 2 3 6 179 32 f - 44 a f + 21 a 180(D16) [k = - -----------------------] 181 2 2 182 256 a 183(C17) step3(2); 184 5 2 5 185 a COS(5 t) + (36 a f - 24 a ) COS(3 t) 186(D17) y (t) = ---------------------------------------- + a SIN(t) + b COS(t) 187 2 1024 2 2 188(C18) step4(2); 189 2 5 190 36 a f - 23 a 191(D18) [[a = 0, b = - ---------------]] 192 2 2 1024 193(C19) step5(2); 194 5 2 5 195 a COS(5 t) + (36 a f - 24 a ) COS(3 t) 196(D19) y (t) = ---------------------------------------- 197 2 1024 198 199 2 5 200 (36 a f - 23 a ) COS(t) 201 - ------------------------ 202 1024 203 204 205The process can be automated 206(C20) duffing(3); 207 3 3 208 a COS(3 t) a COS(t) 209y (t) = ----------- - --------- 210 1 32 32 211 212 5 2 5 2 213 a COS(5 t) 9 a f COS(3 t) 3 a COS(3 t) 9 a f COS(t) 214y (t) = ----------- + --------------- - ------------- - ------------- 215 2 1024 256 128 256 216 217 5 218 23 a COS(t) 219 + ------------ 220 1024 221 222 7 4 7 2 223 a COS(7 t) 13 a f COS(5 t) 3 a COS(5 t) 81 a f COS(3 t) 224y (t) = ----------- + ---------------- - ------------- + ---------------- 225 3 32768 6144 2048 2048 226 227 4 7 2 4 228 423 a f COS(3 t) 297 a COS(3 t) 81 a f COS(t) 1217 a f COS(t) 229 - ----------------- + --------------- - -------------- + ---------------- 230 8192 16384 2048 24576 231 232 7 233 547 a COS(t) 234 - ------------- 235 32768 236 237 3 3 3 2 6 9 2 2 3 6 238 e (128 f - 236 a f + 221 a f - 81 a ) e (32 f - 44 a f + 21 a ) 239w= - ------------------------------------------ - ---------------------------- 240 3 2 241 2048 a 256 a 242 243 3 244 e (4 f - 3 a ) 245 - -------------- + 1 246 8 a 247 248 249The patch applied is 250 251--- duffing.mac.orig 2004-02-14 11:42:26.389030400 +1100 252+++ duffing.mac 2004-02-14 11:27:35.948640000 +1100 253@@ -40,7 +40,8 @@ 254 :expand(trigreduce(expand(ev(eq[i],makelist([e[j],f[j]],j,1,i-1), 255 diff))))$ 256 step2(i):=(f[i]:solve(coeff(temp1,cos(t)),k[i]),temp1:ev(temp1,f[i]))$ 257-step3(i):=temp1:ev(ode2(temp1,y[i](t),t),%k1:a[i],%k2:b[i])$ 258+step3(i):=(temp1:ode2(temp1,y[i](t),t), 259+ temp1:subst(a[i],%k1,temp1),temp1:subst(b[i],%k2,temp1))$ 260 step4(i):=(temp2:rhs(temp1),temp2:diff(temp2,t), 261 temp2:solve([ev(rhs(temp1),t:0),ev(temp2,t:0)],[a[i],b[i]]))$ 262 step5(i):=e[i]:ev(temp1,temp2)$ 263 264 265 266Local Variables: *** 267mode: Text *** 268End: ***