Dobbeltpendul og kaos med maple restart: with(plots): with(plottools): g:=9.82; #Her kan g \303\246ndres, hvis n\303\270dvendigt! Denne maplesheet er et produkt af to timers kedsomhed og derfor ikke super gennemarbejdet. Men jeg tror at det er meget l\303\246rerigt at lege lidt med et dobbeltpendul og udforske kaos n\303\246rmere. Det hele er oprindeligt lavet for eget forn\303\270jelse, men jeg har tilf\303\270jet nogle ord hvis andre var interesseret i at smugkigge.

Herned bliver funktionenerne DoublePendul() og PendulAnimation() defineret til at l\303\270se og animere problemet. Tryk enter og luk denne sektion! #----Definition af differentialligninger---- ligning1:= proc(theta1,theta2,m1,m2,L1,L2) diff(theta1,t,t) = ( -g*(2*m1+m2)*sin(theta1) - g*m2*sin(theta1-2*theta2) - 2*sin(theta1-theta2)*m2*(diff(theta2,t)^2*L2+diff(theta1,t)^2*L1*cos(theta1-theta2)) )/ ( L1*(2*m1+m2-m2*cos(2*theta1-2*theta2)) ) end proc: ligning2:= proc(theta1,theta2,m1,m2,L1,L2) diff(theta2,t,t) = ( 2*sin(theta1-theta2)*(diff(theta1,t)^2*L1*(m1+m2)+g*(m1+m2)*cos(theta1)+diff(theta2,t)^2*L2*m2*cos(theta1-theta2)) )/ ( L2*(2*m1+m2-m2*cos(2*theta1-2*theta2))) end proc: printf("Ligninger der skal l\303\270ses:"); ligning1(theta1(t),theta2(t),m1,m2,L1,L2); ligning2(theta1(t),theta2(t),m1,m2,L1,L2); printf(" "); printf(" "); #----Ligningsl\303\270ser---- DobbeltPendul:= proc(theta1_0,theta2_0,d_theta1_0,d_theta2_0, m1, m2, L1, L2)::listprocedure;local dsys,n_frames; dsys := {ligning1(theta1(t),theta2(t),m1,m2,L1,L2),ligning2(theta1(t),theta2(t),m1,m2,L1,L2),theta1(0)=theta1_0,theta2(0)=theta2_0, D(theta1)(0)=d_theta1_0, D(theta2)(0)=d_theta2_0}: return [op(dsolve(dsys,numeric,output=listprocedure)),masse1 = m1,masse2 = m2, Laengde1 = L1, Laengde2 = L2]: end proc: #----Animation---- PendulAnimation:= proc(InputData,t_end) #--------------------- m1:= subs(InputData,masse1); m2:= subs(InputData,masse2); L1:= subs(InputData,Laengde1); L2:= subs(InputData,Laengde2); n_frames:= t_end*10; m1_radius:=m1*0.1: m2_radius:=m2*0.1: theta_1:=subs(InputData,theta1(t)); theta_2:=subs(InputData,theta2(t)); #--------------------- #----Prototyper af figur---- Circle := (x,y,r)->circle([x,y], r, color=blue): Line := (x1,y1,x2,y2)->line([x1,y1], [x2,y2], color=red,thickness=2): baggrund:=PLOT(rectangle([-L1-L2-1,0],[L1+L2+1,-L1-L2-1],color=white)): #----Tegning af pendul---- m1_centrum:=n->display(Circle(L1*sin(theta_1(0.1*n)),-L1*cos(theta_1(0.1*n)),m1_radius)): m2_centrum:=n->display(Circle(L1*sin(theta_1(0.1*n))+L2*sin(theta_2(0.1*n)),-L1*cos(theta_1(0.1*n))-L2*cos(theta_2(0.1*n)),m2_radius)): m1_line:=n->display(Line(0,0,L1*sin(theta_1(0.1*n)),-L1*cos(theta_1(0.1*n)) )): m2_line:=n->display(Line(L1*sin(theta_1(0.1*n)),-L1*cos(theta_1(0.1*n)),L1*sin(theta_1(0.1*n))+L2*sin(theta_2(0.1*n)),-L1*cos(theta_1(0.1*n))-L2*cos(theta_2(0.1*n)) )): p:=n->display(m1_centrum(n),m2_centrum(n),m1_line(n),m2_line(n)): #Animation pr n'te frame A := array(1..2): A[1]:=animate(p,[n],n=0..n_frames,frames=n_frames,scaling=constrained,background=baggrund): A[2]:=odeplot(InputData,[theta1(t),theta2(t)],t=0..t_end,frames=n_frames,title='Faseplot'): display(A); end proc:

For at l\303\270se ligningerne skal man v\303\246lge nogle startsbetingelser til vinklerne LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEnJiM5NTI7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYkLUkjbW5HRiQ2JFEiMUYnRjVGNS8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnLUkobWZlbmNlZEdGJDYkLUYjNiQtRi82JVEidEYnL0YzUSV0cnVlRicvRjZRJ2l0YWxpY0YnRjVGNUY1 og LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEnJiM5NTI7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYkLUkjbW5HRiQ2JFEiMkYnRjVGNS8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnLUkobWZlbmNlZEdGJDYkLUYjNiQtRi82JVEidEYnL0YzUSV0cnVlRicvRjZRJ2l0YWxpY0YnRjVGNUY1 (se figur 1.4, s. 14 i bogen) samt masser og l\303\246ngder. For at l\303\270se ligningerne numerisk benyttes syntaksen DoublePendulLUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkobWZlbmNlZEdGJDYmLUYjNjgtSSVtc3ViR0YkNiUtSSNtaUdGJDYnUScmIzk1MjtGJy8lJWJvbGRHUSV0cnVlRicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSVib2xkRicvJStmb250d2VpZ2h0R0Y/LUYjNiYtSSNtbkdGJDYmUSIxRidGN0Y9RkBGN0Y9RkAvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1JI21vR0YkNi9RIixGJ0Y3Rj1GQC8lJmZlbmNlR0Y8LyUqc2VwYXJhdG9yR0Y5LyUpc3RyZXRjaHlHRjwvJSpzeW1tZXRyaWNHRjwvJShsYXJnZW9wR0Y8LyUubW92YWJsZWxpbWl0c0dGPC8lJ2FjY2VudEdGPC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYxNiVGMy1GIzYmLUZFNiZRIjJGJ0Y3Rj1GQEY3Rj1GQEZIRkstRjE2JS1JJm1vdmVyR0YkNiUtRiM2JkYzRjdGPUZALUYjNiYtRkw2L1EiLkYnRjdGPUZARk8vRlJGPEZTRlVGV0ZZRmVuRmduL0Zbb0ZpbkY3Rj1GQC8lJ2FjY2VudEdGOUZCRkhGSy1GMTYlLUZnbzYlLUYjNiYtRjQ2J1EoJnRoZXRhO0YnRjdGOkY9RkBGN0Y9RkBGW3BGYnBGX29GSEZLLUZMNi9RIn5GJ0Y3Rj1GQEZPRmBwRlNGVUZXRllGZW5GZ25GYXBGXXFGXXFGXXEtRjE2JS1GNDYnUSJtRidGNy9GO0Y5L0Y+USxib2xkLWl0YWxpY0YnRkBGQkZIRkstRjE2JUZicUZfb0ZIRkstRjE2JS1GNDYnUSJMRidGN0ZlcUZmcUZARkJGSEZLLUYxNiVGXHJGX29GSEY3Rj1GQEY3Rj1GQEY3Rj1GQA==, hvor \316\270'erne er til tiden 0, og som output kommer der en liste med l\303\270sningerne. Lad os starte med et eksempel: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzY+LUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEnJiM5NTI7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYkLUkjbW5HRiQ2JFEiMUYnRjVGNS8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnLUkobWZlbmNlZEdGJDYkLUYjNiQtRjs2JFEiMEYnRjVGNUY1LUkjbW9HRiQ2LVEiPUYnRjUvJSZmZW5jZUdGNC8lKnNlcGFyYXRvckdGNC8lKXN0cmV0Y2h5R0Y0LyUqc3ltbWV0cmljR0Y0LyUobGFyZ2VvcEdGNC8lLm1vdmFibGVsaW1pdHNHRjQvJSdhY2NlbnRHRjQvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0Znbi1JJm1mcmFjR0YkNigtRi82JVEnJiM5NjA7RidGMkY1LUYjNiQtRjs2JFEiMkYnRjVGNS8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGam8vJSliZXZlbGxlZEdGNC1GSjYtUSIsRidGNUZNL0ZQUSV0cnVlRidGUUZTRlVGV0ZZL0ZmblEmMC4wZW1GJy9GaW5RLDAuMzMzMzMzM2VtRictRko2LVEifkYnRjVGTUZPRlFGU0ZVRldGWUZkcC9GaW5GZXBGaHAtRiw2JUYuRmBvRj5GQUZJRl1vRl9wRmhwRmhwRmhwLUYsNiUtSSZtb3ZlckdGJDYlLUYjNiRGLkY1LUYjNiQtRko2L1EiLkYnLyUlYm9sZEdGY3AvRjZRJWJvbGRGJy8lK2ZvbnR3ZWlnaHRHRl1yRk1GT0ZRRlNGVUZXRllGZHBGW3FGNS8lJ2FjY2VudEdGY3BGOEY+RkFGSUZGRl9wRmhwRmhwRmhwLUYsNiUtRmFxNiUtRiM2JC1GLzYlUSgmdGhldGE7RidGMkY1RjVGZXFGYHJGYG9GPkZBRklGRkY1 LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzZALUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEibUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUkjbW5HRiQ2JFEiMUYnL0Y2USdub3JtYWxGJ0Y+LyUvc3Vic2NyaXB0c2hpZnRHUSIwRictSSNtb0dGJDYtUSI9RidGPi8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGSS8lKXN0cmV0Y2h5R0ZJLyUqc3ltbWV0cmljR0ZJLyUobGFyZ2VvcEdGSS8lLm1vdmFibGVsaW1pdHNHRkkvJSdhY2NlbnRHRkkvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZYLUY7NiRRIjVGJ0Y+LUZENi1RIixGJ0Y+RkcvRktGNEZMRk5GUEZSRlQvRldRJjAuMGVtRicvRlpRLDAuMzMzMzMzM2VtRictRkQ2LVEifkYnRj5GR0ZKRkxGTkZQRlJGVEZcby9GWkZdb0Zgb0Zgb0Zgb0Zgby1GLDYlRi4tRiM2JC1GOzYkUSIyRidGPkY+RkBGQy1GOzYkUSI0RidGPkZobkZgb0Zgb0Zgb0Zgby1GLDYlLUYvNiVRIkxGJ0YyRjVGOEZARkNGW3BGaG5GYG9GYG9GYG9GYG9GYG8tRiw2JUZgcEZmb0ZARkNGZW5GPg== losning:= DobbeltPendul(Pi/2,Pi,0,0,5,4,4,5); hvor listen med l\303\270sningen er gemt i losning. For at animere bev\303\246gelsen skal man skrive PendulAnimation(losning, t_end), hvor t_end er l\303\246ngden af animationen (i sekunder). F.eks. (kan tage lidt tid at beregne): PendulAnimation(losning, 15); Klik p\303\245 plottet og der kommer en STOP/PLAY-menu op i toppen (alternativ: h\303\270jreklik-->Animation-->Play). For at hive en af funktionerne ud af listen med l\303\270sninger, kan subs() benyttes. F.eks. hvis vi er interesseret i at arbejde med LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEnJiM5NTI7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYkLUkjbW5HRiQ2JFEiMUYnRjVGNS8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnLUkobWZlbmNlZEdGJDYkLUYjNiQtRi82JVEidEYnL0YzUSV0cnVlRicvRjZRJ2l0YWxpY0YnRjVGNUY1: theta1:=subs(losning, theta1(t)); og denne kan benyttes som en normal funktion, f.eks.: theta1(4); plot(theta1(t),t=0..10); #reset af variable losning:='losning': theta1 :='theta1': Der f\303\270lger nogle eksempler p\303\245 anvendelser.

Her er et simpel eksempel med to n\303\246sten identiske systemer, hvor det eneste forskel er at i det ene system er pendulets LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2J1EnJiM5NTI7RicvJSVib2xkR1EldHJ1ZUYnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1ElYm9sZEYnLyUrZm9udHdlaWdodEdGOi1GIzYmLUkjbW5HRiQ2JlEiMUYnRjJGOEY7RjJGOEY7LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvRjlRJ25vcm1hbEYn startet 0.02 radianer (ca. 1 grad) anderledes end det andet system. Normalt (as in "line\303\246re systemer") ville s\303\245dan en forskel v\303\246re negligibelt. losning1 :=DobbeltPendul(Pi/2, Pi/4, 0, 0, 7,5,7,4): #1 losning2 :=DobbeltPendul(Pi/2+0.02, Pi/4, 0, 0, 7,5,7,4): #2 Og her ses animationerne (hiv i plottet s\303\245 det passer til sk\303\246rmst\303\270rrelsen): #Kan tage lidt tid at beregne Grid:= array(1..2,1..1): Grid[1,1]:=PendulAnimation(losning1,30): Grid[2,1]:=PendulAnimation(losning2,30): display(Grid); optil ca. 15 s f\303\270lger systemerne (n\303\246sten identisk) hinanden, hvor de derefter \303\246ndrer (meget groft) "opf\303\270rsel"! dette kan ses tydeligere ved at se p\303\245 LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEnJiM5NTI7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYkLUkjbW5HRiQ2JFEiMkYnRjVGNS8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnLUkobWZlbmNlZEdGJDYkLUYjNiQtRi82JVEidEYnL0YzUSV0cnVlRicvRjZRJ2l0YWxpY0YnRjVGNUY1 for de to systemer: theta_2_1:=subs(losning1,theta2(t)): #1 theta_2_2:=subs(losning2,theta2(t)): #2 plot([theta_2_1(t),theta_2_2(t)],t=0..30); som det ses er de identiske til ca. t = 15 sek. Et dobbelpendul er derfor meget f\303\270lsom overfor startsbetingelserne! #reset af variable theta_2_2:='theta_2_2': theta_2_1:='theta_2_1': losning1:='losning1': losning2:='losning2':

N\303\245r man taylorudvikler ligningerne og l\303\270ser dem, f\303\245r man p\303\246ne harmoniske l\303\270sninger for sm\303\245 udsving. Her kan vi unders\303\270ge dette. Vi har her et pendul i hvile, der bliver skubbet igang. SmaaUdsving:= DobbeltPendul(0,0,0.1,0.1,8,10,10,4): PendulAnimation(SmaaUdsving,15); Bev\303\246gelsen er som forventet. Lad os plotte vinklerne for funktion af tiden: theta1:=subs(SmaaUdsving,theta1(t)): theta2:=subs(SmaaUdsving,theta2(t)): plot([theta1(t),theta2(t)],t=0..20); Som forventet. Det kunne v\303\246re sk\303\246gt at sammenligne vinkelfrekvensen med de lineariserede l\303\270sninger, dog ikke sk\303\246g nok til at jeg gider det nu. #reset af variable SmaaUdsving:='SmaaUdsving': theta1:='theta1': theta2:='theta2':

There better be! #Koordinater x1:=t->L1*sin(theta1(t)); y1:=t->-L1*cos(theta1(t)); x2:=t->L1*sin(theta1(t))+L2*sin(theta2(t)); y2:=t->-L1*cos(theta1(t))-L2*cos(theta2(t)); #Kinetisk og potential energi V:=t->m1*g*y1(t) + m2*g*y2(t); T:=t->1/2*m1*(D[1](x1)(t)^2+D[1](y1)(t)^2) + 1/2*m2*(D[1](x2)(t)^2+D[1](y2)(t)^2); #Et tilf\303\246ldigt system L1:=7: L2:=6: m1:=6: m2:=7: losning3:= DobbeltPendul(0,Pi/3, 2,0,m1,m2,L1,L2): #PendulAnimation(losning3,15); og s\303\245 plot af energien over tid theta1:=subs(losning3,theta1(t)): theta2:=subs(losning3,theta2(t)): odeplot(losning3,[t, T(t)+V(t) ],t=0..200,axes=boxed,view=[0..200,170..180]); Det ser ud som om det totale energi af systemet vokser (meget langsomt) og ikke er bevaret! Men ingen grund til at tro at alt vores teori er forkert og en nobelpris er i sigte! Dette skyldes at der sker en masse afrundinger n\303\245r ligningerne l\303\270ses numerisk og over langtid kommer der derfor sm\303\245 fejl. N\303\245r nu systemet er s\303\245 f\303\270lsom overfor startsbetingelser, kan det t\303\246nkes at systemet ogs\303\245 er meget f\303\270lsom overfor s\303\245danne sm\303\245 afrundingsfejl. Efter et stykke tid, vil vores l\303\270sning formentlig v\303\246re helt galt p\303\245 den! Derfor er det i praksis umuligt at forudsige bev\303\246gelsen af et dobbeltpendul, selvom det er et s\303\245 simpelt system. S\303\245 hvor sv\303\246rt mon det er at forudsige vejret? #reset af variable x1:='x1': y1:='y1': x2:='x2': y2:='y2': V:='V': T:='T': L1:='L1': L2:='L2': m1:='m1': m2:='m2': theta1:='theta1': theta2:='theta2':

Hvis du alligevel ikke har meget andet at fortage dig her i livet, leg l\303\270s med dobbeltpendul selv! P1:=DobbeltPendul(0,Pi/2,0,0,10,100,50,7): PendulAnimation(P1,50);