Model wahadła matematycznego

Łukasz Błaszczyk
Metody komputerowe w równaniach różniczkowych, ćwiczenie 3

Rozważymy dwa modele wahadła matematycznego - model pojedynczego wahadła oraz model wahadła podwójnego. Znajdziemy równanie ruchu, które opisuje zachowanie wahadła w obu przypadkach. Szczególnie interesuje nas równanie wahadła podwójnego, którego rozwiązanie może dać wyniki dość chaotyczne, tak jak na rysunku poniżej.

g=9.81; m1=1; m2=1; l1=5; l2=3; u10=(π/3); u20=(π/2);
solution=NDSolve[{
(m1+m2)l1 v1'[t]+m2 l2 v2'[t] Cos[u2[t]-u1[t]]==m2 l2 ((v2[t])^2)Sin[u2[t]-u1[t]]-(m1+m2)g Sin[u1[t]],
l2 v2'[t]+l1 v1'[t] Cos[u2[t]-u1[t]]==-l1 ((v1[t])^2)Sin[u2[t]-u1[t]]-g Sin[u2[t]],
v1[t]==u1'[t],
v2[t]==u2'[t],
u1[0]==u10,u2[0]==u20,v1[0]==0,v2[0]==0},{u1[t],u2[t],v1[t],v2[t]},{t,0,100}];
rr1[t_]=Evaluate[{l1 Sin[u1[t]],-l1 Cos[u1[t]]}/.solution][[1]];
rr2[t_]=rr1[t]+Evaluate[{l2 Sin[u2[t]],-l2 Cos[u2[t]]}/.solution][[1]];
Manipulate[
Show[
Graphics[{LightGray,EdgeForm[],Rectangle[{-1,0},{1,1}]}],
Graphics[{Black,Thick,Line[{{-1,0},{1,0}}]}],
Graphics[{Black,Dashed,Line[{{0,0},{0,-9}}]}],
Graphics[{Orange,Thick,Line[{{0,0},rr1[t],rr2[t]}]}],
ParametricPlot[rr2[s],{s,0,t},PlotStyle->Directive[Green,Thick]],
Graphics[{PointSize[0.03],Orange,Point[rr1[t]],Point[rr2[t]]}],
PlotRange-> {{-8,8},{-10,1}}
],{t,0.01,100},SaveDefinitions->True]

Tym razem model (z fizycznego punktu widzenia) jest dużo prostszy niż w przypadku równania wiszącego łańcucha. Skorzystamy jedynie z II Zasady Dynamiki Newtona.

Wahadło matematyczne

Rozważmy najprostszy przypadek, tzn. pojedynczy punk o masie m zawieszony na nieważkiej, nierozciągliwej nici o długości L. Niech θ oznacza kąt między osią pionową, a nicią.

Show[
Graphics[{LightGray,EdgeForm[],Rectangle[{-1,0},{1,1}]}],
Graphics[{Black,Thick,Line[{{-1,0},{1,0}}]}],
Graphics[{Black,Dashed,Line[{{0,0},{0,-8}}]}],
Graphics[{Orange,Thick,Line[{{0,0},{1.5*2.5,-1.5*4.330127018922194`}}]}],
Graphics[{PointSize[0.03],Red,Point[{1.5*2.5,-1.5*4.330127018922194`}]}],
Graphics[{Text[Style["θ",Large],{0.6,-2},FormatType->StandardForm]}],
Graphics[{Text[Style["L",Large],{2.8,-3.2},FormatType->StandardForm]}],
Graphics[{Text[Style["m",Large],{2.8,-6.7},FormatType->StandardForm]}],
ParametricPlot[{3 Sin[t],-3 Cos[t]},{t,0,(π/6)},PlotStyle->Black],
PlotRange-> {{-8,8},{-10,1}}
]

Zauważmy, że ruch wahadła wystarczy opisać za pomocą kąta θ zmieniającego się w czasie. Wówczas ruch punktu (w układzie kartezjańskim) opisuje funkcja wektorowa

Clear[g,L,m]

R[t_]=L{Sin[u[t]],-Cos[u[t]]};

wsp=0.4;
P={1.5*2.5,-1.5*4.330127018922194`};
Q={-1.5*4.330127018922194`,-1.5*2.5};
PQx = -P[[1]]/Q[[1]];
Show[
Graphics[{LightGray,EdgeForm[],Rectangle[{-1,0},{1,1}]}],
Graphics[{Black,Thick,Line[{{-1,0},{1,0}}]}],
Graphics[{Black,Dashed,Line[{{0,0},{0,-8}}]}],
Graphics[{Orange,Thick,Line[{{0,0},{1.5*2.5,-1.5*4.330127018922194`}}]}],
Graphics[{Blue,Thick,Arrowheads[0.06],Arrow[{P,P-(wsp/PQx)P}]}],
Graphics[{Blue,Thick,Arrowheads[0.06],Arrow[{P,P+wsp(Q+(P/PQx))}]}],
Graphics[{Black,Dashed,Line[{P-(wsp/PQx)P,P+wsp Q,P+wsp(Q+(P/PQx))}]}],
Graphics[{Darker[Green],Thick,Arrowheads[0.06],Arrow[{P,P+wsp Q}]}],
Graphics[{Text[Style["-(T/L)R[t]",Large,Blue],{4.8,-3.6},FormatType->StandardForm]}],
Graphics[{Text[Style["F",Large,Darker[Green]],{2.5,-6.7},FormatType->StandardForm]}],
Graphics[{Text[Style["{0,-mg}",Large,Blue],{6.1,-9.7},FormatType->StandardForm]}],
Graphics[{PointSize[0.03],Red,Point[{1.5*2.5,-1.5*4.330127018922194`}]}],
PlotRange-> {{-8,8},{-12,1}}
]

Jakie siły działają w tym układzie? Mamy siłę ciężkości działającą w dół (równą {0,-mg}) oraz naprężenie nici o wartości T i kierunku -R[t]/L. Wypadkowa siła jest więc równa

F[t_]=-(T/L)R[t]+{0,-m g};

Korzystając z II Zasady Dynamiki Newtona otrzymamy więc

eq=m R''[t]-F[t]

i powyższy wektor musi być zerowy.

Zauważmy, że mamy dwa równania i... dwie niewiadome. Nie znamy również naprężenia nici. Spróbujmy jednak wyeliminować ten parametr z równania - zauważmy, że wystarczy pomnożyć pierwsze równanie przez Cos[u[t]], a drugie pomnożyć przez Sin[u[t]] i dodać oba równania do siebie.

TrigReduce[eq[[1]] Cos[u[t]]+eq[[2]] Sin[u[t]]]

Widzimy, że masa jest tutaj nieistotnym parametrem. Takie równanie możemy rozwiązać numerycznie (musimy podać warunek brzegowy na funkcję u[t] oraz na jej pierwszą pochodną w zerze), jednak najpierw zbadajmy portret fazowy tego równania.

g=9.81; L=8;

StreamPlot[{v,-(g/L)Sin[u]},{u,-2π,2π},{v,-2π,2π}]

Można wyciągnąć z tego wnioski na temat tego jak będzie wyglądałorozwiązanie w zależności od początkowego kąta oraz początkowej prędkości jaką nadamy wahadłu.

solution=NDSolve[{g Sin[u[t]]+L u''[t]==0,u[0]==(π/3),u'[0]==0},u[t],{t,0,20}];

r[t_]=Evaluate[{L Sin[u[t]],-L Cos[u[t]]}/.solution][[1]];

Manipulate[
Show[
Graphics[{LightGray,EdgeForm[],Rectangle[{-1,0},{1,1}]}],
Graphics[{Black,Thick,Line[{{-1,0},{1,0}}]}],
Graphics[{Black,Dashed,Line[{{0,0},{0,-8}}]}],
Graphics[{Orange,Thick,Line[{{0,0},r[t]}]}],
ParametricPlot[r[s],{s,0,t},PlotStyle->Directive[Green,Thick]],
Graphics[{PointSize[0.03],Red,Point[r[t]]}],
PlotRange-> {{-8,8},{-10,1}}
],{t,0.01,20},SaveDefinitions->True]

Wahadło podwójne

W przypadku podwójnego wahadła sytuacja jest analogiczna. Tym razem mamy dwa punkty o masach m1 i m2 na niciach o długościach L1 i L2. Wprowadzamy dwie funkcje kątowe θ1 i θ2, które w pełni opisują ruch wahadła.

Show[
Graphics[{LightGray,EdgeForm[],Rectangle[{-1,0},{1,1}]}],
Graphics[{Black,Thick,Line[{{-1,0},{1,0}}]}],
Graphics[{Black,Dashed,Line[{{0,0},{0,-8}}]}],
Graphics[{Black,Dashed,Line[{{2.5,-4.330127018922194`},{2.5,-8}}]}],
Graphics[{Orange,Thick,Line[{{0,0},{2.5,-4.330127018922194`},{5.098076211353316`,-5.8301270189221945`}}]}],
Graphics[{PointSize[0.03],Red,Point[{2.5,-4.330127018922194`}],Point[{5.098076211353316`,-5.8301270189221945`}]}],
Graphics[{Text[Style["(Subscript[θ, 1])",Large],{0.7,-2},FormatType->StandardForm]}],
Graphics[{Text[Style["(Subscript[θ, 2])",Large],{3.3,-5.3},FormatType->StandardForm]}],
Graphics[{Text[Style["(Subscript[L, 1])",Large],{2.2,-2.2},FormatType->StandardForm]}],
Graphics[{Text[Style["(Subscript[L, 2])",Large],{4.2,-4.4},FormatType->StandardForm]}],
Graphics[{Text[Style["(Subscript[m, 1])",Large],{1.7,-4.2},FormatType->StandardForm]}],
Graphics[{Text[Style["(Subscript[m, 2])",Large],{6.2,-6.0},FormatType->StandardForm]}],
ParametricPlot[{3 Sin[t],-3 Cos[t]},{t,0,(π/6)},PlotStyle->Black],
ParametricPlot[{2.5,-4.330127018922194`}+{2 Sin[t],-2 Cos[t]},{t,0,(π/3)},PlotStyle->Black],
PlotRange-> {{-8,8},{-10,1}}
]

Clear[g,L1,L2,m1,m2]

R1[t_]=L1 {Sin[u1[t]],-Cos[u1[t]]};
R2[t_]=R1[t]+L2{Sin[u2[t]],-Cos[u2[t]]};

wsp=2;
wsp1=0.5;
wsp2=0.8;
P={2.5,-4.330127018922194`};
Q={-4.330127018922194`,-2.5};
R={5.098076211353316`,-5.8301270189221945`};
S={R[[2]]-P[[2]],P[[1]]-R[[1]]};
PQx=-P[[1]]/Q[[1]];
RSx=-(R[[1]]-P[[1]])/S[[1]];

vec1=(wsp2/RSx)(R-P);
Show[
Graphics[{LightGray,EdgeForm[],Rectangle[{-1,0},{1,1}]}],
Graphics[{Black,Thick,Line[{{-1,0},{1,0}}]}],
Graphics[{Black,Dashed,Line[{{0,0},{0,-8}}]}],
Graphics[{Black,Dashed,Line[{P,{2.5,-8}}]}],
Graphics[{Orange,Thick,Line[{{0,0},P,R}]}],
Graphics[{Blue,Thick,Arrowheads[0.03],Arrow[{P,P-1.7vec1[[1]](P/(P[[1]]))}]}],
Graphics[{Blue,Thick,Arrowheads[0.03],Arrow[{P,P-0.7vec1[[1]] (Q/(Q[[1]]))+1.7vec1[[1]](P/(P[[1]]))-(wsp2/RSx)(R-P)}]}],
Graphics[{Blue,Thick,Arrowheads[0.03],Arrow[{P,P+(wsp2/RSx)(R-P)}]}],
Graphics[{Green,Dashed,Arrowheads[0.02],Arrow[{P,P-0.7vec1[[1]] (Q/(Q[[1]]))-(wsp2/RSx)(R-P)}]}],
Graphics[{Black,Dashed,Line[{P-1.7vec1[[1]](P/(P[[1]])),P-0.7vec1[[1]] (Q/(Q[[1]]))-(wsp2/RSx)(R-P), P-0.7vec1[[1]] (Q/(Q[[1]]))+1.7vec1[[1]](P/(P[[1]]))-(wsp2/RSx)(R-P)}]}],
Graphics[{Black,Dashed,Line[{ P-0.7vec1[[1]] (Q/(Q[[1]]))-(wsp2/RSx)(R-P),P-0.7vec1[[1]] (Q/(Q[[1]])),P+(wsp2/RSx)(R-P)}]}],
Graphics[{Darker[Green],Thick,Arrowheads[0.03],Arrow[{P,P-0.7vec1[[1]] (Q/(Q[[1]]))}]}],

Graphics[{Magenta,Thick,Arrowheads[0.03],Arrow[{R,R-(wsp2/RSx)(R-P)}]}],
Graphics[{Magenta,Thick,Arrowheads[0.03],Arrow[{R,R+wsp2(S+((R-P)/RSx))}]}],
Graphics[{Black,Dashed,Line[{R-(wsp2/RSx)(R-P),R+wsp2 S,R+wsp2(S+((R-P)/RSx))}]}],
Graphics[{Brown,Thick,Arrowheads[0.03],Arrow[{R,R+wsp2 S}]}],

Graphics[{Text[Style["-((Subscript[T, 1])/(Subscript[L, 1]))(Subscript[R, 1])[t]",Medium,Blue],{3.1,-3.2},FormatType->StandardForm]}],
Graphics[{Text[Style["(Subscript[F, 1])",Medium,Darker[Green]],{1.6,-4.4},FormatType->StandardForm]}],
Graphics[{Text[Style["{0,-(Subscript[m, 1])g}",Medium,Blue],{1.7,-7.8},FormatType->StandardForm]}],
Graphics[{Text[Style["((Subscript[T, 2])/(Subscript[L, 1]))((Subscript[R, 2])[t]-(Subscript[R, 1])[t])",Medium,Blue],{4.9,-4.4},FormatType->StandardForm]}],
Graphics[{Text[Style["((Subscript[T, 2])/(Subscript[L, 1]))((Subscript[R, 1])[t]-(Subscript[R, 2])[t])",Medium,Magenta],{6.2,-5.1},FormatType->StandardForm]}],
Graphics[{Text[Style["(Subscript[F, 2])",Medium,Brown],{3.7,-7.6},FormatType->StandardForm]}],
Graphics[{Text[Style["{0,-(Subscript[m, 2])g}",Medium,Magenta],{6.0,-8.8},FormatType->StandardForm]}],
Graphics[{PointSize[0.03],Red,Point[P],Point[R]}],

PlotRange-> {{-8,8},{-10,1}},ImageSize->Large
]

W podwójnym układzie siły musimy rozpatrywać osobno dla każdego punktu. W przypadku punktu o masie m1 mamy siłę ciężkości oraz dwa naprężenia - jedno ciągnie punkt w górę (T1), natomiast drugie w dół (T2). Na punkt o masie m2 działa siła ciężkości oraz jedno naprężenie w górę, równe naprężeniu w dół dla punktu o masie m1 (czyli T2).

F1[t_]=-(T1/L1)R1[t]+(T2/L2)(R2[t]-R1[t])+{0,-m1 g};
F2[t_]=-(T2/L2)(R2[t]-R1[t])-m2{0,g};

Ponownie korzystamy z II Zasady Dynamiki Newtona, tym razem jednak dla obu punktów.

eq1=m1 R1''[t]-F1[t]

eq2=m2 R2''[t]-F2[t]

i powyższe wektory muszą być zerowe.

Podobnie jak wcześniej mamy nadmiarowe zmienne (czyli funkcje T1 i T2). Spróbujmy się ich jakoś pozbyć. Zobaczmy, że mnożąc pierwsze równanie przez Cos[u1[t]] i drugie przez Sin[u1[t]] i dodając je do siebie pozbędziemy się funkcji T1. Jednocześnie przy u1’’[t] zredukują się wyrazy trygonometryczne (z jedynki trygonometrycznej). Aby pozbyć się funkcji T2 pomnóżmy trzecie równanie przez Cos[u1[t]], czwartek przez Sin[u1[t]] i dodajmy wszystko do siebie. W ten sposób otrzymamy jedno równanie (drugiego rzędu) wiążące funkcje u1 i u2.

EQ1=FullSimplify[TrigReduce[eq1[[1]] Cos[u1[t]]+eq1[[2]] Sin[u1[t]]+eq2[[1]] Cos[u1[t]]+eq2[[2]]Sin[u1[t]]]/.{u1'[t]-> v1[t],u1''[t]-> v1'[t],u2'[t]-> v2[t],u2''[t]-> v2'[t]}]

Potrzeba jednak również drugiego równania. Aby je otrzymać, wystarczy trzecie równanie pomnożyć przez Cos[u2[t]], czwarte przez Sin[u2[t]] i dodać te dwa równania do siebie. Znikną wtedy również funkcje trygonometryczne przy drugich pochodnych i równanie się uprości.

EQ2=FullSimplify[TrigReduce[eq2[[1]]Cos[u2[t]]+eq2[[2]]Sin[u2[t]]]/.{u1'[t]-> v1[t],u1''[t]-> v1'[t],u2'[t]-> v2[t],u2''[t]-> v2'[t],m2-> 1}]

Mając te dwa równania możemy już je rozwiązać. Tym razem badanie portretu fazowego jest raczej nierealne.

g=9.81; L1=3; L2=2; m1=1; m2=2;

solution=NDSolve[{EQ1==0,EQ2==0,v1[t]==u1'[t],v2[t]==u2'[t],
u1[0]==(π/6),u2[0]==(π/2),v1[0]==0,v2[0]==0},{u1[t],u2[t],v1[t],v2[t]},{t,0,20}];

Plot[Evaluate[{u1[t],u2[t]}/.solution],{t,0,20},PlotLegends->{"u1[t]","u2[t]"}]

Widać z powyższego wykresu, że ruch takiego wahadła nie jest okresowy - więcej szczegółów widać na animacji wahadła.

r1[t_]=Evaluate[{L1 Sin[u1[t]],-L1 Cos[u1[t]]}/.solution][[1]];
r2[t_]=r1[t]+Evaluate[{L2 Sin[u2[t]],-L2 Cos[u2[t]]}/.solution][[1]];

Manipulate[
Show[
Graphics[{LightGray,EdgeForm[],Rectangle[{-1,0},{1,1}]}],
Graphics[{Black,Thick,Line[{{-1,0},{1,0}}]}],
Graphics[{Black,Dashed,Line[{{0,0},{0,-9}}]}],
Graphics[{Orange,Thick,Line[{{0,0},r1[t],r2[t]}]}],
ParametricPlot[r2[s],{s,0,t},PlotStyle->Directive[Green,Thick]],
Graphics[{PointSize[0.03],Orange,Point[r1[t]],Point[r2[t]]}],
PlotRange-> {{-6,6},{-7,1}}
],{t,0.01,20},SaveDefinitions->True]

Dobierając różne parametry długości (i różne warunki początkowe) można pokazać chaotyczne zachowanie wahadła.