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

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,4*Pi/3},AxesLabel->{"x [m]","y [m]"}]

Sin[Pi/4]

Cos[Pi/4]

sinw=Sqrt[(1-1/Sqrt[2])/2];

cosw=Simplify[Sqrt[1-sinw^2]];

N[sinw-Sin[Pi/8]]

N[cosw-Cos[Pi/8]]

(* xO i yO odpowiadają u= Pi/8 *)

xO=r*(Pi/8-sinw);

yO=-r*(1-cosw);

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

xP=x1/.{u->4*Pi/3}

yP=y1/.{u->4*Pi/3}

pkty={{xO,yO},{xP,yP}};

N[pkty]

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

Show[plot0,plot1,PlotRange->All]

linia prosta: y=a*x+b
Teraz parametr u ma sens położenia punktu na osi x.

a=(yP-yO)/(xP-xO);

b=yO-a*xO;

x2=u;

y2=a*u+b;

plot2=ParametricPlot[{x2,y2},{u,xO,xP},PlotStyle->{Green,Dashed}]

Show[plot0,plot1,plot2,PlotRange->All]

Parabola nr 1
Parametr u ma sens położenia punktu na osi x.

sa012=Simplify[Solve[{a2*xO^2+a1*xO+a0==yO,a2*xP^2+a1*xP+a0==yP,-a1/(2*a2)==-2/5},{a0,a1,a2}]];

a0=a0/.sa012[[1]];

a1=a1/.sa012[[1]];

a2=a2/.sa012[[1]]

x3=u;

y3=a0+a1*u+a2*u^2;

plot3=ParametricPlot[{x3,y3},{u,0,xP},PlotStyle->{Brown,Dotted,Thick}]

Show[plot0,plot1,plot2,plot3,PlotRange->All]

Parabola nr 2
Parametr u ma sens położenia punktu na osi x.

d=(yO-yP)/(xO-xP)^2;

x4=u;

y4=d*(u-xP)^2+yP;

plot4=ParametricPlot[{x4,y4},{u,0,xP},PlotStyle->{Gray}]

Show[plot0,plot1,plot2,plot3,plot4,PlotRange->All]

Parabola nr 3
Parametr u ma sens położenia punktu na osi x.

snewa012=Simplify[Solve[{b2*xO^2+b1*xO+b0==yO,b2*xP^2+b1*xP+b0==yP,b2*3^2+b1*3+b0==-3},{b0,b1,b2}]];

b0=b0/.snewa012[[1]];

b1=b1/.snewa012[[1]];

b2=b2/.snewa012[[1]];

x5=u;

y5=b0+b1*u+b2*u^2;

plot5=ParametricPlot[{x5,y5},{u,0,xP},PlotStyle->{Black,Dashed,Thick}]

Show[plot0,plot1,plot2,plot3,plot4,plot5,PlotRange->All]

Teraz modelujemy ruch koralika po wybranych krzywych. 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 *)

Zaczynamy od brachistochrony.

ClearAll[x1,y1];

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

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

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

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

(* Energia kinetyczna punktowej masy *)

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

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

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

(* Funkcja Lagrange'a układu *)

L=FullSimplify[Ekin-Epot];

Etot=FullSimplify[Ekin+Epot];

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

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

d3=D[d1,t];

s1=Simplify[Solve[d3==d2,f1''[t]]];

s1[[1]];

eq1=s1[[1]]/. Rule->Equal;

(* Warunki początkowe: koralik znajduje się w punkcie O(xO,yO), a jego początkowa prędkość wynosi zero. Oznacza to, że f1[0]=Pi/8 oraz f1'[0]=0    *)

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

(* maksymalny czas śledzenia ruchu koralika *)

tmax=3; (* s *)

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

Flatten[numerycznie1];

(* Wyciągam funkcję f1[t] z numerycznego rozwiązania *)

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

(* Teraz mogę używać funkcji f1[t] tak, jak innych funkcji *)

f1[1]

f1'[0.5]

Plot[f1[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","f1(t) [rad]"}]

Plot[f1'[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","f1'(t) [rad/s]"}]

Plot[Evaluate[Ekin],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Ekin(t) [J]"}]

Plot[Evaluate[Epot],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Epot(t) [J]"}]

Plot[Evaluate[Etot],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Etot(t) [J]"}]

end1=FindRoot[f1[t]-4*Pi/3,{t,1}]

t1=t/.end1

x1[t_]:=Evaluate[(r*(f1[t]-Sin[f1[t]]))/.Flatten[numerycznie1]];

y1[t_]:=Evaluate[(-r*(1-Cos[f1[t]]))/.Flatten[numerycznie1]];

ParametricPlot[{x1[t],y1[t]},{t,0,t1},PlotRange->All]

Animate[Show[plot1,ListPlot[{{x1[t],y1[t]}},PlotStyle->{PointSize[0.03],Orange}]],{t,0,t1}]

Teraz rozważamy ruch po prostej

ClearAll[x2,y2];

x2[t]=f2[t];

y2[t]=a*f2[t]+b;

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

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

(* Energia kinetyczna punktowej masy *)

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

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

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

(* Funkcja Lagrange'a układu *)

L=FullSimplify[Ekin-Epot];

Etot=FullSimplify[Ekin+Epot];

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

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

d3=D[d1,t];

s2=Simplify[Solve[d3==d2,f2''[t]]];

s2[[1]];

eq2=s2[[1]]/. Rule->Equal;

(* Warunki początkowe: koralik znajduje się w punkcie O(xO,yO), a jego początkowa prędkość wynosi zero. Oznacza to, że f2[0]=xO oraz f2'[0]=0    *)

wp2={f2[0]==xO,f2'[0]==0};

(* maksymalny czas śledzenia ruchu koralika *)

tmax=3; (* s *)

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

Flatten[numerycznie2];

(* Wyciągam funkcję f2[t[ z numerycznego rozwiązania *)

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

(* Teraz mogę używać funkcji f2[t] tak, jak innych funkcji *)

f2[1]

f2'[0.5]

Plot[f2[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","f2(t) [m]"}]

Plot[f2'[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","f2'(t) [rad/s]"}]

Plot[Evaluate[Ekin],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Ekin(t) [J]"}]

Plot[Evaluate[Epot],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Epot(t) [J]"}]

Plot[Evaluate[Etot],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Etot(t) [J]"}]

end2=FindRoot[f2[t]-xP,{t,1}]

t2=t/.end2

x2[t_]:=Evaluate[(f2[t])/.Flatten[numerycznie2]];

y2[t_]:=Evaluate[(a*f2[t]+b)/.Flatten[numerycznie2]];

ParametricPlot[{x2[t],y2[t]},{t,0,t2},PlotRange->All]

Animate[Show[plot2,ListPlot[{{x2[t],y2[t]}},PlotStyle->{PointSize[0.03],Orange}]],{t,0,t2}]

Teraz rozważamy ruch po pierwszej paraboli

ClearAll[x3,y3];

x3[t]=f3[t];

y3[t]=a0+a1*f3[t]+a2*(f3[t])^2;

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

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

(* Energia kinetyczna punktowej masy *)

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

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

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

(* Funkcja Lagrange'a układu *)

L=Simplify[Ekin-Epot];

Etot=Simplify[Ekin+Epot];

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

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

d3=D[d1,t];

s3=Simplify[Solve[d3==d2,f3''[t]]];

eq3=s3[[1]]/. Rule->Equal;

(* Warunki początkowe: koralik znajduje się w punkcie O(xO,yO), a jego początkowa prędkość wynosi zero. Oznacza to, że f3[0]=xO oraz f'[0]=0    *)

wp3={f3[0]==xO,f3'[0]==0};

(* maksymalny czas śledzenia ruchu koralika *)

tmax=5; (* s *)

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

Flatten[numerycznie3];

(* Wyciągam funkcję f3[t[ z numerycznego rozwiązania *)

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

(* Teraz mogę używać funkcji f3[t] tak, jak innych funkcji *)

f3[1]

f3'[0.5]

Plot[f3[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","f3(t) [m]"}]

Plot[f3'[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","f3'(t) [m/s]"}]

Plot[Evaluate[Ekin],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Ekin(t) [J]"}]

Plot[Evaluate[Epot],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Epot(t) [J]"}]

Plot[Evaluate[Etot],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Etot(t) [J]"}]

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

t3=t/.end3

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

y3[t_]:=Evaluate[(a0+a1*f3[t]+a2*(f3[t])^2)/.Flatten[numerycznie3]];

ParametricPlot[{x3[t],y3[t]},{t,0,t3},PlotRange->All]

Animate[Show[plot3,ListPlot[{{x3[t],y3[t]}},PlotStyle->{PointSize[0.03],Orange}]],{t,0,t3}]

Teraz rozważamy ruch po drugiej paraboli

ClearAll[x4,y4];

x4[t]=f4[t];

y4[t]=d*(f4[t]-xP)^2+yP;

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

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

(* Energia kinetyczna punktowej masy *)

Ekin=Simplify[(1/2)*m*(x4'[t]^2+y4'[t]^2)];

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

Epot=  Simplify[m*g*y4[t]];

(* Funkcja Lagrange'a układu *)

L=Simplify[Ekin-Epot];

Etot=Simplify[Ekin+Epot];

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

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

d3=D[d1,t];

s4=Simplify[Solve[d3==d2,f4''[t]]];

eq4=s4[[1]]/. Rule->Equal;

(* Warunki początkowe: koralik znajduje się w punkcie O(xO,yO), a jego początkowa prędkość wynosi zero. Oznacza to, że f4[0]=xO oraz f4'[0]=0    *)

wp4={f4[0]==xO,f4'[0]==0};

(* maksymalny czas śledzenia ruchu koralika *)

tmax=3; (* s *)

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

Flatten[numerycznie4];

(* Wyciągam funkcję f4[t[ z numerycznego rozwiązania *)

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

(* Teraz mogę używać funkcji f4[t] tak, jak innych funkcji *)

f4[1]

f4'[0.5]

Plot[f4[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","f4(t) [m]"}]

Plot[f4'[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","f4'(t) [m/s]"}]

Plot[Evaluate[Ekin],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Ekin(t) [J]"}]

Plot[Evaluate[Epot],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Epot(t) [J]"}]

Plot[Evaluate[Etot],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Etot(t) [J]"}]

end4=FindRoot[f4[t]-xP,{t,1}]

t4=t/.end4

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

y4[t_]:=Evaluate[(d*(f4[t]-xP)^2+yP)/.Flatten[numerycznie4]];

ParametricPlot[{x4[t],y4[t]},{t,0,t4},PlotRange->All]

Animate[Show[plot4,ListPlot[{{x4[t],y4[t]}},PlotStyle->{PointSize[0.03],Orange}]],{t,0,t4}]

Teraz rozważamy ruch po trzeciej paraboli

ClearAll[x5,y5];

x5[t]=f5[t];

y5[t]=b2*(f5[t])^2+b1*f5[t]+b0;

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

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

(* Energia kinetyczna punktowej masy *)

Ekin=Simplify[(1/2)*m*(x5'[t]^2+y5'[t]^2)];

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

Epot=  Simplify[m*g*y5[t]];

(* Funkcja Lagrange'a układu *)

L=Simplify[Ekin-Epot];

Etot=Simplify[Ekin+Epot];

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

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

d3=D[d1,t];

s5=Simplify[Solve[d3==d2,f5''[t]]];

eq5=s5[[1]]/. Rule->Equal;

(* Warunki początkowe: koralik znajduje się w punkcie O(xO,yO), a jego początkowa prędkość wynosi zero. Oznacza to, że f5[0]=xO oraz f5'[0]=0    *)

wp5={f5[0]==xO,f5'[0]==0};

(* maksymalny czas śledzenia ruchu koralika *)

tmax=5; (* s *)

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

Flatten[numerycznie5];

(* Wyciągam funkcję f5[t[ z numerycznego rozwiązania *)

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

(* Teraz mogę używać funkcji f5[t] tak, jak innych funkcji *)

f5[1]

f5'[0.5]

Plot[f5[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","f5(t) [m]"}]

Plot[f5'[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","f5'(t) [m/s]"}]

Plot[Evaluate[Ekin],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Ekin(t) [J]"}]

Plot[Evaluate[Epot],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Epot(t) [J]"}]

Plot[Evaluate[Etot],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Etot(t) [J]"}]

end5=FindRoot[f5[t]-xP,{t,1}]

t5=t/.end5

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

y5[t_]:=Evaluate[(b2*(f5[t])^2+b1*f5[t]+b0)/.Flatten[numerycznie5]];

ParametricPlot[{x5[t],y5[t]},{t,0,t5},PlotRange->All]

Animate[Show[plot5,ListPlot[{{x5[t],y5[t]}},PlotStyle->{PointSize[0.03],Orange}]],{t,0,t5}]

x11[t_]=If[t≤t1,Evaluate[x1[t]],xP];

y11[t_]=If[t≤t1,Evaluate[y1[t]],yP];

x21[t_]=If[t≤t2,Evaluate[x2[t]],xP];

y21[t_]=If[t≤t2,Evaluate[y2[t]],yP];

x31[t_]=If[t≤t3,Evaluate[x3[t]],xP];

y31[t_]=If[t≤t3,Evaluate[y3[t]],yP];

x41[t_]=If[t≤t4,Evaluate[x4[t]],xP];

y41[t_]=If[t≤t4,Evaluate[y4[t]],yP];

x51[t_]=If[t≤t5,Evaluate[x5[t]],xP];

y51[t_]=If[t≤t5,Evaluate[y5[t]],yP];

Animate[Show[plot1,plot2,plot3,plot4,plot5,ListPlot[{{x11[t],y11[t]}},PlotStyle->{PointSize[0.03],Orange}],
ListPlot[{{x21[t],y21[t]}},PlotStyle->{PointSize[0.03],Red}],
ListPlot[{{x31[t],y31[t]}},PlotStyle->{PointSize[0.03],Green}],
ListPlot[{{x41[t],y41[t]}},PlotStyle->{PointSize[0.03],Brown}],
ListPlot[{{x51[t],y51[t]}},PlotStyle->{PointSize[0.03],Pink}],PlotRange->All
],{t,0,t3}]

OSTATECZNE WYNIKI: WYGRYWA BRACHISTOCHRONA

{t1,t2,t3,t4,t5}