W tym notebooku przekonamy się, że czas ruchu po brachistochronie z położenia początkowego (wybranego powyżej punktu najniższego) do jej najniższego punktu nie zależy od wyboru punktu początkowego !

Polecam także stronę:
https://kierul.wordpress.com/2014/05/16/johann-bernoulli-i-brachistochrona-krzywa-najszybszego-spadku-1697/



(obrazek z powyżej polecanej strony)

brachistochrona (fragment cykloidy)
W równaniach parametrycznych występuje wielkość u.
Dla brachistochrony u ma sens kąta.

x1=r*(u-Sin[u]);

y1=-r*(1-Cos[u]);

r=1;plot1=ParametricPlot[{x1,y1},{u,0,Pi},AxesLabel->{"x [m]","y [m]"}]

(* xO1 i yO1 odpowiadają u= Pi/8 *)

xO1=x1/.{u-> Pi/8};

yO1=y1/.{u-> Pi/8};

(* xO2 i yO2 odpowiadają u= Pi/5 *)

xO2=x1/.{u-> Pi/5};

yO2=y1/.{u-> Pi/5};

(* xO3 i yO3 odpowiadają u= Pi/3 *)

xO3=x1/.{u-> Pi/3};

yO3=y1/.{u-> Pi/3};

(* xO4 i yO4 odpowiadają u= 2*Pi/3 *)

xO4=x1/.{u-> 2*Pi/3};

yO4=y1/.{u-> 2*Pi/3};

(* xO4 i yO4 odpowiadają u= 2*Pi/3 *)

xO4=x1/.{u-> 2*Pi/3};

yO4=y1/.{u-> 2*Pi/3};

(* xO5 i yO5 odpowiadają u= 5*Pi/6 *)

xO5=x1/.{u-> 5*Pi/6};

yO5=y1/.{u-> 5*Pi/6};

(* xP i yP odpowiadają u=Pi *)

xP=x1/.{u->Pi};

yP=y1/.{u->Pi};

pkty={{xO1,yO1},{xO2,yO2},{xO3,yO3},{xO4,yO4},{xO5,yO5},{xP,yP}};

N[pkty]

plot0=ListPlot[pkty,PlotStyle->{Red,PointSize[Large]},PlotRange->All]

Show[plot0,plot1,PlotRange->All]

ClearAll[x1,y1];

Teraz modelujemy ruch koralika po wybranej brachistochronie, wybierając różne położenia początkowe, ale zakładając, że początkowa prędkość koralika wynosi zero.

Naszym celem jest wyznaczenie, jak zmienia się w funkcji czasu parametr określający położenie koralika od czasu.

(* Wartość przyspieszenia ziemskiego *)

g=981/100; (* m/s^2 *)

(* Ruch koralika bez tarcia nie zależy od jego masy, ale ponieważ liczymy jego energię kinetyczną i potencjalną, przyjmijmy jakąś wartość m *)

m=1/100; (* kg *)

ClearAll[x,y];

x[t]=r*(f[t]-Sin[f[t]]);

y[t]=-r*(1-Cos[f[t]]);

x'[t]=D[x[t],t];

y'[t]=D[y[t],t];

(* Energia kinetyczna punktowej masy *)

Ekin=FullSimplify[(1/2)*m*(x'[t]^2+y'[t]^2)];

(* Energia potencjalna zależy wyłącznie od y1 *)

Epot=  FullSimplify[m*g*y[t]];

(* Funkcja Lagrange'a układu *)

L=FullSimplify[Ekin-Epot];

Etot=FullSimplify[Ekin+Epot];

d1=D[L,f'[t]];

d2=D[L,f[t]];

d3=D[d1,t];

s=Simplify[Solve[d3==d2,f''[t]]];

s[[1]];

eq=s[[1]]/. Rule->Equal

eq1=eq/.f->f1

eq2=eq/.f->f2

eq3=eq/.f->f3

eq4=eq/.f->f4

eq5=eq/.f->f5

(* Warunki początkowe: koralik znajduje się w punkcie Oi(xOi,yOi), a jego początkowa prędkość wynosi zero.   *)

wp1={f1[0]==Pi/8,f1'[0]==0};

wp2={f2[0]==Pi/5,f2'[0]==0};

wp3={f3[0]==Pi/3,f3'[0]==0};

wp4={f4[0]==2*Pi/3,f4'[0]==0};

wp5={f5[0]==5*Pi/6,f5'[0]==0};

(* maksymalny czas śledzenia ruchu koralika *)

tmax=50; (* s *)

numerycznie1=NDSolve[Flatten[{eq1,wp1}],f1[t],{t,0,tmax},WorkingPrecision->24];

f1[t_]:=Evaluate[f1[t]/.Flatten[numerycznie1]]

x1[t_]=r*(f1[t]-Sin[f1[t]]);

y1[t_]=-r*(1-Cos[f1[t]]);

numerycznie2=NDSolve[Flatten[{eq2,wp2}],f2[t],{t,0,tmax},WorkingPrecision->24];

f2[t_]:=Evaluate[f2[t]/.Flatten[numerycznie2]]

x2[t_]=r*(f2[t]-Sin[f2[t]]);

y2[t_]=-r*(1-Cos[f2[t]]);

numerycznie3=NDSolve[Flatten[{eq3,wp3}],f3[t],{t,0,tmax},WorkingPrecision->24];

f3[t_]:=Evaluate[f3[t]/.Flatten[numerycznie3]]

x3[t_]=r*(f3[t]-Sin[f3[t]]);

y3[t_]=-r*(1-Cos[f3[t]]);

numerycznie4=NDSolve[Flatten[{eq4,wp4}],f4[t],{t,0,tmax},WorkingPrecision->24];

f4[t_]:=Evaluate[f4[t]/.Flatten[numerycznie4]]

x4[t_]=r*(f4[t]-Sin[f4[t]]);

y4[t_]=-r*(1-Cos[f4[t]]);

numerycznie5=NDSolve[Flatten[{eq5,wp5}],f5[t],{t,0,tmax},WorkingPrecision->24];

f5[t_]:=Evaluate[f5[t]/.Flatten[numerycznie5]]

x5[t_]=r*(f5[t]-Sin[f5[t]]);

y5[t_]=-r*(1-Cos[f5[t]]);

Sprawdzimy czasy ruchu koralików puszczonych z różnych wysokości.
Najpierw wyznaczymy całkowite czasy ruchów: t1, t2 i t3.

end1=FindRoot[f1[t]-Pi,{t,1}]//Quiet

t1=t/.end1

end2=FindRoot[f2[t]-Pi,{t,1}]//Quiet

t2=t/.end2

end3=FindRoot[f3[t]-Pi,{t,1}]//Quiet

t3=t/.end3

end4=FindRoot[f4[t]-Pi,{t,1}]//Quiet

t4=t/.end4

end5=FindRoot[f5[t]-Pi,{t,1}]//Quiet

t5=t/.end5

Animate[Show[plot1,ParametricPlot[Evaluate[{x1[t],y1[t]}],{t,0.9999*s,s},PlotStyle->{PointSize[0.03],Orange}]/.Line->Point,
ParametricPlot[Evaluate[{x2[t],y2[t]}],{t,0.9999*s,s},PlotStyle->{PointSize[0.03],Red}]/.Line->Point,
ParametricPlot[Evaluate[{x3[t],y3[t]}],{t,0.9999*s,s},PlotStyle->{PointSize[0.03],Green}]/.Line->Point,
ParametricPlot[Evaluate[{x4[t],y4[t]}],{t,0.9999*s,s},PlotStyle->{PointSize[0.03],Brown}]/.Line->Point,
ParametricPlot[Evaluate[{x5[t],y5[t]}],{t,0.9999*s,s},PlotStyle->{PointSize[0.03],Blue}]/.Line->Point
],{s,0,t1}]

aaa=Table[Show[plot1,ParametricPlot[Evaluate[{x1[t],y1[t]}],{t,0.9999*s,s},PlotStyle->{PointSize[0.03],Orange}]/.Line->Point,
ParametricPlot[Evaluate[{x2[t],y2[t]}],{t,0.9999*s,s},PlotStyle->{PointSize[0.03],Red}]/.Line->Point,
ParametricPlot[Evaluate[{x3[t],y3[t]}],{t,0.9999*s,s},PlotStyle->{PointSize[0.03],Green}]/.Line->Point,
ParametricPlot[Evaluate[{x4[t],y4[t]}],{t,0.9999*s,s},PlotStyle->{PointSize[0.03],Brown}]/.Line->Point,
ParametricPlot[Evaluate[{x5[t],y5[t]}],{t,0.9999*s,s},PlotStyle->{PointSize[0.03],Blue}]/.Line->Point
],{s,0.001,t1,0.2}];

Export[
 NotebookDirectory[] <> "brachistochrona_2.gif", aaa]