N-násobné kyvadlo

nKyv=5;  (*počet hmotných bodů a závěsů*)

Pojmenování veličin


x[1],y[1],...x[5],y[5]  ... polohy hmotných bodů
ϕ[1][t] ... ϕ[5][t]  ... Zobecněné souřadnice problému
l[1],m[1],...l[5],m[5]  ... parametry hmotných bodů a závěsů

1. Výpočet lagranžiánu

{x[0],y[0]} = {0,0}; (* x[0],y[0] jsou pomocné veličiny *)
Do[
{x[k] ,y[k]} = {x[k-1] ,y[k-1]} + l[k]{Sin[ϕ[k][t]],-Cos[ϕ[k][t]]}
,{k,1,nKyv}]

{x[5],y[5]}

eT = Sum[m[k]/2( D[x[k],t]^2+D[y[k],t]^2),{k,1,nKyv}];
eV = Sum[g m[k] y[k],{k,1,nKyv}];
eL = eT-eV;

2. Pohybové rovnice

LR2[L_,q_]:= D[ D[L,q'[t]] ,t]-D[L,q[t]]==0

rovnice = Table[LR2[ eL, ϕ[k]],{k,1,nKyv}];

3. Parametry pro numerické řešení

numParam=(*{m1->1,m2->1,g->1,l1->1,l2->1,ϕ0->Pi/2}*)
Join[{ϕ0->Pi/2,g->1},Table[{l[k]->1,m[k]->1},{k,1,nKyv}]//Flatten]

pocPodm=Table[{ϕ[k][0]==ϕ0,ϕ[k]'[0]==0},{k,1,nKyv}]/.numParam//Flatten

4. Explicitní pohybové rovnice jsou nepoužitelné

(podobně neumí NDSolve řešit soustavu s nevyjádřenými nejvyššími derivacemi)

(*
explicitniRovnice=Solve[ rovnice,Table[ϕ[k]''[t],{k,1,nKyv}] ][[1]]/.Rule->Equal;

explicitniRovnice/.numParam/.(Table[{ϕ[k][t]->Pi/2,ϕ[k]'[0]->0},{k,1,nKyv}]//Flatten)
*)

(*
tMax=44;
Join[ explicitniRovnice , pocPodm ]/.numParam;
numReseni = NDSolve[%,{ϕ1[t],ϕ2[t]},{t,0,tMax}][[1]]
*)

5. Numerické řešení oklikou

matSoustavy=D[ rovnice[[All,1]], {Table[ (ϕ[k])''[t],{k,1,nKyv}]}]/.numParam;
pravaStrana= -rovnice[[All,1]]/. Table[ (ϕ[k])''[t]->0,{k,1,nKyv}]/.numParam;

ϕ2[ϕ0:{__Real},ϕ1:{__Real}]:=Module[{},
subst = Table[{ (ϕ[k])'[t]->ϕ1[[k]],ϕ[k][t]->ϕ0[[k]]},{k,1,nKyv}]//Flatten;
LinearSolve[ matSoustavy/.subst, pravaStrana/.subst ]
]

(* kontrola *)
ϕ2[Table[Pi/2,nKyv]//N,Table[0,nKyv]//N]

tMax=84;
pocUhly=Table[ϕ0,nKyv]/.numParam//N;
pocRychlosti=Table[0,nKyv]//N;

qFce=NDSolve[ D[q[t],t,t]==ϕ2[q[t],D[q[t],t]] &&
q'[0]==pocRychlosti && q[0]==pocUhly,q[t],{t,0,tMax}][[1,1,2,0]];

casyT=Range[0,tMax,0.05];
nCasu=Length@casyT;
numUhly = Table[ qFce[t],{t,casyT}];
casyRules[iT_Integer]:=(Table[ϕ[k][t]->numUhly[[iT,k]],{k,0,nKyv}]);
numPolohy=Table[ Table[{x[k],y[k]},{k,0,nKyv}]/.numParam/.casyRules[iT],{iT,1,nCasu}];

Table[
gZavazi = Point[numPolohy[[it]]];
gZavesy = Line[numPolohy[[it]]];
Rasterize@Graphics[{Blue,PointSize[0.02],gZavazi,Red,gZavesy},AspectRatio->Automatic,PlotRange->{{-nKyv,nKyv},{-nKyv,nKyv/3}}],
{it,1,nCasu}]; (*Table + ListAnimate uslozi animaci do *.nb *)

ListAnimate[%,25]

Table[
it0=Max[1,it-25];
gStiny = Table[Line[Table[numPolohy[[it2,k+1]],{it2,it0,it}]],{k,1,nKyv}];
gZavazi = Point[numPolohy[[it]]];
gZavesy = Line[numPolohy[[it]]];
Rasterize@Graphics[{Lighter@Lighter@Gray,Thickness[0.005],gStiny,Red,gZavesy,Blue,PointSize[0.02],gZavazi},AspectRatio->Automatic,PlotRange->{{-nKyv,nKyv},{-nKyv,nKyv/3}}],
{it,1,nCasu}];
ListAnimate[%,25]