1 /*
2 * $Id: demo3.i,v 1.1 2005-09-18 22:05:55 dhmunro Exp $
3 * Chaotic pendulum demo
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 func demo3(time_limit, kick=, damp_time=)
12 /* DOCUMENT demo3
13 Solve Lagrange's equations for a famous chaotic pendulum (as on
14 exhibit at the San Francisco Exploratorium). Run a movie of the
15 result. By default, the movie runs for 60 seconds, but if you
16 supply an argument to the demo3 function, it will run for that
17 many seconds.
18
19 The kick= keyword may be used to adjsut the initial amplitude;
20 kick= 1.2 is the default.
21
22 You may also wish to supply a damp_time= keyword to see the effect
23 of an ad hoc damping term. damp_time=100 is nice.
24 */
25 {
26 require, "rkutta.i";
27 a= 2.0;
28 b= 3.0;
29 c= 1.5;
30 gravity= 1.0;
31 widget_setup, a, b, c, 1.00, 0.67, gravity;
32
33 dt= 0.05*sqrt(min(a,b,c)/gravity);
34
35 if (is_void(damp_time)) damp_time= 1.e30;
36 if (is_void(kick)) kick= 1.2;
37
38 /* roll dice to get initial velocities */
39 thetadot= kick;
40 gammadot= 0.1*random();
41 qdotq= [[thetadot,0.,0.,gammadot], [0.,0.,0.,0.]];
42
43 window, wait=1;
44 s= 1.1*(b+c);
45 limits, -s,s,-s,s;
46 display, qdotq; /* without this, won't get axis numbers */
47
48 require, "movie.i";
49 movie, draw_frame, time_limit, 0.02;
50 }
51
draw_frame(i)52 func draw_frame(i)
53 {
54 display, qdotq;
55 qdotq= rkdumb(lagrange, qdotq,0., dt, 1, nosave=1);
56 return 1; /* just go until time limit expires */
57 }
58
widget_setup(a,b,c,m_tee,m_bar,gravity)59 func widget_setup(a, b, c, m_tee, m_bar, gravity)
60 {
61 /* set up the values used by lagrange */
62 extern Tdiag, Ta, Tb, Tc, Ba, Bb, Bg, Dg;
63
64 B= 0.5*m_bar*c^2;
65 C= 0.5*B*c;
66 A= m_tee*(2.*a^3+b^3)/3. + m_bar*c*(2.*a^2+b^2) + C;
67 Ba= B*a;
68 Bb= B*b;
69
70 Dg= (0.5*m_tee*b^2 + m_bar*b*c)*gravity;
71 Bg= B*gravity;
72
73 Tdiag= [[A,0.,0.,0.], [0.,C,0.,0.], [0.,0.,C,0.], [0.,0.,0.,C]];
74 Ta= [[0.,Ba,0.,0.], [Ba,0.,0.,0.], [0.,0.,0.,0.], [0.,0.,0.,0.]];
75 Tb= [[0.,0.,Ba,0.], [0.,0.,0.,0.], [Ba,0.,0.,0.], [0.,0.,0.,0.]];
76 Tc= [[0.,0.,0.,Bb], [0.,0.,0.,0.], [0.,0.,0.,0.], [Bb,0.,0.,0.]];
77 }
78
lagrange(qdotq,time)79 func lagrange(qdotq, time/*unused*/)
80 {
81 /* compute the time derivatives of [qdot,q]
82 Lagrange's equations d(dL/dqdot)/dt - dL/dq = 0
83 turn out to be T*qdotdot+L=0, with T and L computed here. */
84 theta= qdotq(1,2);
85 alpha= qdotq(2,2);
86 beta= qdotq(3,2);
87 gamma= qdotq(4,2);
88 qdot= qdotq(,1);
89 t2= qdot(1)^2;
90 a2= qdot(2)^2;
91 b2= qdot(3)^2;
92 g2= qdot(4)^2;
93
94 amt= alpha-theta;
95 bpt= beta+theta;
96 gmt= gamma-theta;
97
98 T= Tdiag + Ta*sin(amt) + Tb*sin(bpt) - Tc*cos(gmt);
99 camt= cos(amt);
100 cbpt= cos(bpt);
101 sgmt= sin(gmt);
102 L= [Dg*sin(theta) + Ba*(a2*camt+b2*cbpt)+Bb*g2*sgmt,
103 Bg*sin(alpha) - Ba*t2*camt,
104 Bg*sin(beta) + Ba*t2*cbpt,
105 Bg*sin(gamma) - Bb*t2*sgmt];
106 qdotdot= LUsolve(T, -L);
107
108 qdotdot(,1)-= qdot/damp_time;
109
110 return [qdotdot, qdot];
111 }
112
display(qdotq)113 func display(qdotq)
114 {
115 /* draw the widget given its current state */
116 theta= qdotq(1,2);
117 alpha= qdotq(2,2);
118 beta= qdotq(3,2);
119 gamma= qdotq(4,2);
120
121 st= sin(theta);
122 ct= cos(theta);
123 p1x= -a*ct;
124 p1y= -a*st;
125 p2x= -p1x;
126 p2y= -p1y;
127 p3x= b*st;
128 p3y= -b*ct;
129 q1x= p1x + c*sin(alpha);
130 q1y= p1y - c*cos(alpha);
131 q2x= p2x + c*sin(beta);
132 q2y= p2y - c*cos(beta);
133 q3x= p3x + c*sin(gamma);
134 q3y= p3y - c*cos(gamma);
135
136 plg, [p1y,p2y], [p1x,p2x], marks=0,type=1,width=36,color="red";
137 plg, [0.,p3y], [0.,p3x], marks=0,type=1,width=36,color="red";
138 plg, [p1y,q1y], [p1x,q1x], marks=0,type=1,width=24,color="green";
139 plg, [p2y,q2y], [p2x,q2x], marks=0,type=1,width=24,color="green";
140 plg, [p3y,q3y], [p3x,q3x], marks=0,type=1,width=24,color="green";
141 }
142