| blob | 82a515f6c609e40fc64f0276ca5014c9f407de7c |
1 duffing.mac is from the book "Computer Algebra in Applied Mathematics:
2 An introduction to MACSYMA", by Richard H Rand, Pitman (1984).
4 Duffing's equation is
6 x''(T) + x + e*x^3 = e*f*cos(w*T)
8 It is a model of a nonlinear oscillator driven by a sinusiodal forcing
9 function. When e is small and the forcing frequency is close to the
10 natural frequency (of unity) the equation has one or more steady state
11 periodic solutions.
13 This program determines the periodic solutions using Lindstedt's method.
14 The code is very similar to vandpol.mac and mathieu.mac.
16 The original code produces an error. A small patch (below) is required
17 to get the code to run with maxima-5.9.0cvs.
19 The run below, using maxima-5.9.0cvs, reproduce the result on pages
20 127-134 of the book. In the first part, the problem is solved step by
21 step. Then the whole procedure is automated.
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
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
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);
79 Inconsistent equations: (1 2)
80 #0: step4(i=1)
81 -- an error. Quitting. To debug this try DEBUGMODE(TRUE);)
82 (C13)
85 After a small patch
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
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
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
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
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
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
199 2 5
200 (36 a f - 23 a ) COS(t)
201 - ------------------------
202 1024
205 The process can be automated
206 (C20) duffing(3);
207 3 3
208 a COS(3 t) a COS(t)
209 y (t) = ----------- - ---------
210 1 32 32
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)
214 y (t) = ----------- + --------------- - ------------- - -------------
215 2 1024 256 128 256
217 5
218 23 a COS(t)
219 + ------------
220 1024
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)
224 y (t) = ----------- + ---------------- - ------------- + ----------------
225 3 32768 6144 2048 2048
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
232 7
233 547 a COS(t)
234 - -------------
235 32768
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 )
239 w= - ------------------------------------------ - ----------------------------
240 3 2
241 2048 a 256 a
243 3
244 e (4 f - 3 a )
245 - -------------- + 1
246 8 a
249 The patch applied is
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)$
266 Local Variables: ***
267 mode: Text ***
268 End: ***