Autor : Robert Muszyński KCiR PWr Pl  Wszelkie prawa zastrzeżone

Pojedyncze wahadło portret i dywergencja

Model

eqns={ q'[t]==p[t],p'[t]==-g/l Sin[q[t]]-µ/l p[t]} (* µ - współczynnik tarcia *)

(* położenie punktów pośrednich wahadła do wizualizacji *)
pm={l Sin[q[t]],-l Cos[q[t]]};
baza={0,0};
wahadlo={baza,pm};

(* parametry, warunki początkowe, czas symulacji *)
params={l->1,g->1,µ->0};
initial={q[0]==Pi/3,p[0]==0};
tMax=20;
(* rozwiązujemy i wizualizujemy *)
sol=First[NDSolve[{eqns,initial}/.params,{q,p},{t,0,tMax}]];
nFrames=300;
frames=Table[Graphics[{Gray,Thick,Line[wahadlo],Darker[Red],Ball[pm,0.1]}/.params/.sol,PlotRange->{{-1.5,1.5},{-1.5,1.5}},Background->Lighter[White],Axes->True,AxesLabel->{"x","y"}],{t,tMax/nFrames,tMax,tMax/nFrames}];
ListAnimate[frames]

Portret fazowy, pole wektorowe

(* Portret fazowy *)
fC[q_,p_]=p^2-2 g/l Cos[q]/.params;
phase=ContourPlot[fC[q,p],{q,-10,10},{p,-4,4},Contours->{-1.995,-1.9,-1.6,-1,0,1,1.6,1.999,2.7,4,6,8.4,11,14},ContourShading->None,AspectRatio->0.3,ImageSize->Full]

(* Pole wektorowe *)
field=StreamPlot[{p,-g/l Sin[q]-µ/l p}/.params,{q,-10,10},{p,-4,4},StreamPoints->{Automatic,0.4}, StreamColorFunction->"Rainbow",AspectRatio->0.3,ImageSize->Full]

(* Razem *)
Show[field,phase]

(* Orbity na tle portretu fazowego *)
Manipulate[Show[phase,
ParametricPlot[(Evaluate[First[{q[t],p[t]} /.(poss= NDSolve[{eqns,Thread[{q[0],p[0]}==point]}/.params,{q,p},{t,0,T}])]]),{t,0,T},PlotStyle->Red],
wahbase={-9,3};
Graphics[{White,Rectangle[{-10.15,1.8},{-7.85,4.2}],Gray,Thick,Line[wahadlo+{wahbase,wahbase}],Darker[Red],Ball[pm+wahbase,0.1]}/.params/.poss/.t->T,PlotRange->{{-1.5,1.5},{-1.5,1.5}},Background->LightBlue,Axes->True,AxesLabel->{"x","y"}]
],
{{T,0.01},0.01,Tmax},{{Tmax,10},5,40},{{point,{1.05,0}},Locator}, SaveDefinitions->True]

(* Orbity na tle pola wektorowego *)
Manipulate[Show[field,ParametricPlot[Evaluate[First[{q[t],p[t]} /. NDSolve[{eqns,Thread[{q[0],p[0]}==point]}/.params,{q,p},{t,0,T}]]],{t,0,T},PlotStyle->Red]
],{{T,1},0.1,20},{{point,{1.05,0}},Locator}, SaveDefinitions->True]

(* Orbity na tle pola wektorowego - przypadek z tarciem *)fieldfri[fric_]:=StreamPlot[{p,-g/l Sin[q]-µ/l p}/.µ->fric/.params,{q,-10,10},{p,-4,4},StreamPoints->{Automatic,0.4}, StreamMarkers->None,AspectRatio->0.3,ImageSize->Full]
Manipulate[Show[fieldfri[fric],
ParametricPlot[Evaluate[First[{q[t],p[t]} /. (poss=NDSolve[{eqns,Thread[{q[0],p[0]}==point]}/.µ->fric/.params,{q,p},{t,0,T}])]],{t,0,T},PlotStyle->Red]
],
{{T,5},0.1,20},{{fric,0},0,3},{{point,{1.05,0}},Locator}, SaveDefinitions->True]

,,Dywergencja zero”

(* Dywergencja na tle portretu fazowego *)
(*  *)
side[p_,step_:0.1]:=Table[p[[1]]+a(p[[2]]-p[[1]]),{a,0,1,step}]~Join~Table[p[[2]]+a(p[[3]]-p[[2]]),{a,0,1,step}]~Join~Table[p[[3]]+a(p[[1]]-p[[3]]),{a,0,1,step}]
Manipulate[Show[phase,
Graphics[{Red,PointSize[0.002],Point[Evaluate[First[{q[t],p[t]} /. NDSolve[{eqns,Thread[{q[0],p[0]}==#]}/.params,{q,p},{t,0,T}]]&/@side[point1,0.01]/.t->T]]}]],
{{T,0.005},0.005,15},{{point1,{{0,0},{1.05,0},{0,1}}},Locator}, SaveDefinitions->True]

Wahadło Furuty

Wyprowadzenie modelu

Needs["VariationalMethods`"]
Clear[l,a,b,m,θ,ϕ,t,g];
zmienne={θ[t],ϕ[t]}; (* uogólnione współrzędne konfiguracyjne *)
vg= {0,0,-g}; (* wektor grawitacji *)
pm={(a+b Sin[ϕ[t]]) Cos[θ[t]],  (* położenia masy wahadła *)
      (a+b Sin[ϕ[t]])Sin[θ[t]],
       l-b Cos[ϕ[t]]};
vm=D[pm,t]; (* prędkość masy wahadła *)
EK=1/2 m vm. vm; (* Energia kinetyczna *)
EV=-m  pm.vg; (* Energia potencjalna *)
lagranzian=EK-EV;
fureqns=EulerEquations[lagranzian,zmienne,t]

(* położenie punktów pośrednich wahadła do wizualizacji *)
P2=pm/.b->0;
P1=P2/.a->0;
baza={0,0,0};
furuta=Line[{baza,P1,P2,pm}];

Symulacja

(* parametry, warunki początkowe, czas symulacji *)
params={l->1,a->1,b->0.8,m->2,g->9.81};
initial={θ[0]==0,ϕ[0]==0,θ'[0]==0,ϕ'[0]==0};
tMax=20;
(* rozwiązujemy i wizualizujemy *)sol=First[NDSolve[{fureqns,initial}/.params,Head/@zmienne,{t,0,tMax}]];
nFrames=300;
frames=Table[Graphics3D[{Gray,Thick,furuta,Darker[Red],Ball[pm,m/20]}/.params/.sol,PlotRange->{{-2,2},{-2,2},{0,2}},Background->Lighter[White],Boxed->False,Axes->True,AxesLabel->{"x","y","z"}],{t,tMax/nFrames,tMax,tMax/nFrames}];
ListAnimate[frames]

(* trajektoria ruchu wahadła *)
Plot[Evaluate[zmienne/.sol],{t,0,tMax},AxesLabel->{t,zmienne},PlotLegends->zmienne,PlotRange->All]

Podwójne wahadło

Wyprowadzenie modelu

Needs["VariationalMethods`"]
Clear[l1,l2,m1,m2,θ,ϕ,t,g];
zmienne={ϕ1[t],ϕ2[t]}; (* uogólnione współrzędne konfiguracyjne *)
vg= {0,-g}; (* wektor grawitacji *)
(* położenia mas wahadła i wizualizacja wahadła *)
pm1=l1 {Sin[ϕ1[t]],-Cos[ϕ1[t]]};
pm2=pm1+l2 {Sin[ϕ2[t]],-Cos[ϕ2[t]]};
baza={0,0};
pendulum=Line[{baza,pm1,pm2}];
vm1=D[pm1,t]; (* prędkości mas wahadła *)
vm2=D[pm2,t];
(* lagranżian i równania *)
EK=1/2 m1 vm1. vm1 + 1/2 m2 vm2.vm2; (* Energia kinetyczna *)
EV=-m1  pm1.vg - m2 pm2.vg; (* Energia potencjalna *)
lagranzian=EK-EV;
dobeqns=EulerEquations[lagranzian,zmienne,t]

Symulacja

(* parametry, warunki początkowe, czas symulacji *)
params={l1->1/2,l2->1,m1->2,m2->1,g->9.81};
initial={ϕ1[0]==Pi/2,ϕ2[0]==Pi,ϕ1'[0]==0,ϕ2'[0]==0};
tMax=20;
(* rozwiązujemy i wizualizujemy *)sol=First[NDSolve[{dobeqns,initial}/.params,Head/@zmienne,{t,0,tMax}]];
nFrames=500;
frames=Table[Graphics[{Gray,Thick,pendulum,Darker[Blue],Disk[{0,0},.1],Darker[Red],Disk[pm1,m1/20],Disk[pm2,m2/20]}/.params/.sol,PlotRange->{{-2.5,2.5},{-2.5,2.5}},Background->Lighter[White]],{t,0,tMax,tMax/nFrames}];
ListAnimate[frames]

(* trajektoria ruchu wahadła *)
Plot[Evaluate[zmienne/.sol],{t,0,tMax},AxesLabel->{t,zmienne},PlotLegends->zmienne]

(* siła tarcia w pierwszym przegubie - proporcjonalna do prędkości *)
dobeqns[[1,2]]=ϕ1'[t];
dobeqns