Nejprve varování, že býchozí nastavení NDSolve způsobuje nepřesné řešení dif. rovnic s nespojitostmi

NDSolve[{
y''[t]==-y[t],
y[0]==1,
y'[0]==0,
WhenEvent[y[t]==0,y'[t]->-y'[t]]},
y,{t,0,10}]

Plot[y[t]/.%,{t,0,10}]

NDSolve[{
y'[t]==v[t],
v'[t]==-y[t],
y[0]==1,
v[0]==0,
WhenEvent[y[t]==0,v[t]->-v[t]]},
{y[t],v[t]},{t,0,10}]

{y[t],v[t]}/.%[[1]]
Plot[%,{t,0,10}]

Plot[{ Sin[x],Cos[x]},{x,0,10}]

NDSolve[{
y'[t]==v[t],
v'[t]==-y[t],
y[0]==1,
v[0]==0,
WhenEvent[y[t]==0,v[t]->-v[t]]},
{y[t],v[t]},{t,0,10}]

{y[t]^2+v[t]^2-1}/.%[[1]]
Plot[%,{t,0,10}]

reseni=NDSolve[{
y'[t]==v[t],
v'[t]==-y[t],
y[0]==1,
v[0]==0,
WhenEvent[y[t]==0,v[t]->-v[t]]},
{y[t],v[t]},{t,0,10},AccuracyGoal->14,MaxStepSize->0.01][[1]]

energie=y[t]^2+v[t]^2-1/.reseni
Plot[%,{t,0,10}]

casyT = reseni[[1,2,0,3,1]];
energieT = Table[ energie/.t->tk, {tk,casyT}];

ListLinePlot[ Transpose[{casyT,energieT}]]

reseni[[1]]//FullForm;

Rekapitulace

reseni=NDSolve[{
y'[t]==v[t],
v'[t]==-y[t],
y[0]==1,
v[0]==0,
WhenEvent[y[t]==0,v[t]->-v[t]]},
{y[t],v[t]},{t,0,10},AccuracyGoal->14,MaxStepSize->0.01][[1]]

energie=y[t]^2+v[t]^2-1/.reseni
Plot[%,{t,0,10}]

casyT = reseni[[1,2,0,3,1]];
energieT = Table[ energie/.t->tk, {tk,casyT}];

ListLinePlot[ Transpose[{casyT,energieT}]]

Explicitni RK umožní dosáhnot mnohem menší chyby

reseni=NDSolve[{
y'[t]==v[t],
v'[t]==-y[t],
y[0]==1,
v[0]==0,
WhenEvent[y[t]==0,v[t]->-v[t]]},
{y[t],v[t]},{t,0,10},(*StartingStepSize->0.1,*)MaxStepSize->0.01,
 Method->"ExplicitRungeKutta"
][[1]]

energie=y[t]^2+v[t]^2-1/.reseni
Plot[%,{t,0,10}]

casyT = reseni[[1,2,0,3,1]];
energieT = Table[ energie/.t->tk, {tk,casyT}];

ListLinePlot[ Transpose[{casyT,energieT}]]

Explicitni RK + vysoký řád metody

reseni=NDSolve[{
y'[t]==v[t],
v'[t]==-y[t],
y[0]==1,
v[0]==0,WhenEvent[y[t]==0,v[t]->-v[t]]},
{y[t],v[t]},{t,0,10},MaxStepSize->0.1,
Method->{"FixedStep", Method->{"ExplicitRungeKutta",DifferenceOrder->8}}
][[1]]

energie=y[t]^2+v[t]^2-1/.reseni
(*Plot[%,{t,0,10}]*)

casyT = reseni[[1,2,0,3,1]];
energieT = Table[ energie/.t->tk, {tk,casyT}];

ListLinePlot[ Transpose[{casyT,energieT}],PlotLabel->"POcet kroku ="<>ToString[Length[casyT]]]

Problém: disk v homogenním poli uvnitř kružnice

Explicitni RK ale stačí 2.řád, protože volný pád je jen kvadratická funkce času

a = 1
b =1/16;

reseni=NDSolve[{
x'[t]==vx[t],
y'[t]==vy[t],
vx'[t]==0,
vy'[t]==-1,
x[0]==1/10,y[0]==0,vx[0]==0,vy[0]==0,
WhenEvent[   y[t]^2+x[t]^2==(a-b)^2,
{vx[t],vy[t]}->{vx[t]-2x[t] (x[t]vx[t]+y[t]vy[t])/(a-b)^2,vy[t]-2y[t] (x[t]vx[t]+y[t]vy[t])/(a-b)^2}
]
},
{x[t],y[t],vx[t],vy[t]},{t,0,10},MaxStepSize->0.1,
Method->{"FixedStep", Method->{"ExplicitRungeKutta",DifferenceOrder->2}}
][[1]]


ParametricPlot[{x[t],y[t]}/.reseni,{t,0,10},
Epilog->Circle[{0,0},a],PlotRange->{{-a,a},{-a,a}}]


energie=vx[t]^2/2+vy[t]^2/2+y[t]/.reseni;
casyT = reseni[[1,2,0,3,1]];
energieT = Table[ energie/.t->tk, {tk,casyT}];

ListLinePlot[ Transpose[{casyT,energieT}],PlotLabel->"POcet kroku ="<>ToString[Length[casyT]]]

Na Lissajouseových křivkách se naučíme animaci pohybu

xyF[t_]:={Cos[5t],Sin[9t]}
ParametricPlot[xyF[t],{t,0,2Pi}]

Graphics[
{Yellow, (*Circle[{0,0},a],*)Red,Disk[{0,0},b]},
Background->Black,PlotRange->{{-a-b,a+b},{-a-b,a+b}}
]

Graphics[
{Yellow, (*Circle[{0,0},a],*)
RGBColor[0.8,0,0],Disk[{0.05,0.05},b],
Red,Disk[{0,0},b]
},
Background->Black,PlotRange->{{-a-b,a+b},{-a-b,a+b}}
]

Graphics[
{Yellow, (*Circle[{0,0},a],*)
Red,Opacity[0.5],Disk[{0.05,0.05},b],
Opacity[1],Disk[{0,0},b]
},
Background->Black,PlotRange->{{-a-b,a+b},{-a-b,a+b}}
]

a1 ~Plus~ b1

dt=0.05;
stopa = Table[
{Opacity[v],Disk[xyF[3.5+dt*v],b]},{v,0.,1,1/24}];
Graphics[{Red} ~Join~ stopa,
Background->Black,PlotRange->{{-a-b,a+b},{-a-b,a+b}}
]

Animate[
dt=0.05;
stopa = Table[
{Opacity[v],Disk[xyF[T+dt*v],b]},{v,0.,1,1/24}];
Graphics[{Red} ~Join~ stopa,
Background->Black,PlotRange->{{-a-b,a+b},{-a-b,a+b}}
],{T,0,2Pi}]

movie = Table[
dt=0.05;
stopa = Table[
{Opacity[v],Disk[xyF[T+dt*v],b]},{v,0.,1,1/24}];
Graphics[{Red} ~Join~ stopa,
Background->Black,PlotRange->{{-a-b,a+b},{-a-b,a+b}}
],{T,0,2Pi,2Pi/240}];
ListAnimate[movie]

xyF2[t_] := xyF[t+Pi/4]
xyF3[t_] := xyF[t-Pi/4]

movie = Table[
dt=0.05;
stopa = Table[
{Opacity[v],Disk[xyF[T+dt*v],b]},{v,0.,1,1/24}];

stopa2 = Table[
{Opacity[v],Disk[xyF2[T+dt*v],b]},{v,0.,1,1/24}];
stopa3 = Table[
{Opacity[v],Disk[xyF3[T+dt*v],b]},{v,0.,1,1/24}];


Graphics[Join[
{Red}  ,stopa ,
{Blue} , stopa2,
{Darker@Green} , stopa3
],
Background->Black,PlotRange->{{-a-b,a+b},{-a-b,a+b}}
],{T,0,2Pi,2Pi/240}];
ListAnimate[movie]

a = 1
b =1/16;

reseni=NDSolve[{
x'[t]==vx[t],
y'[t]==vy[t],
vx'[t]==0,
vy'[t]==-1,
x[0]==1/10,y[0]==0,vx[0]==0,vy[0]==0,
WhenEvent[   y[t]^2+x[t]^2==(a-b)^2,
{vx[t],vy[t]}->{vx[t]-2x[t] (x[t]vx[t]+y[t]vy[t])/(a-b)^2,vy[t]-2y[t] (x[t]vx[t]+y[t]vy[t])/(a-b)^2}
]
},
{x[t],y[t],vx[t],vy[t]},{t,0,100},MaxStepSize->0.1,
Method->{"FixedStep", Method->{"ExplicitRungeKutta",DifferenceOrder->2}}
][[1]];
xyG[t_Real]:={reseni[[1,2,0]][t],reseni[[2,2,0]][t]}

xyG[94.3]

getF[x0_Real]:=Module[{reseni},
reseni=NDSolve[{
x'[t]==vx[t],
y'[t]==vy[t],
vx'[t]==0,
vy'[t]==-1,
x[0]==x0,y[0]==0,vx[0]==0,vy[0]==0,
WhenEvent[   y[t]^2+x[t]^2==(a-b)^2,
{vx[t],vy[t]}->{vx[t]-2x[t] (x[t]vx[t]+y[t]vy[t])/(a-b)^2,vy[t]-2y[t] (x[t]vx[t]+y[t]vy[t])/(a-b)^2}
]
},
{x[t],y[t],vx[t],vy[t]},{t,0,100},MaxStepSize->0.1,
Method->{"FixedStep", Method->{"ExplicitRungeKutta",DifferenceOrder->2}}
][[1]];
Function[{t},{reseni[[1,2,0]][t],reseni[[2,2,0]][t]}]
]

tmp =getF[0.2]

tmp[3.4]

xyG=getF[0.1002];
xyG2=getF[0.1003];
xyG3=getF[0.1004];

movie = Table[
dt=0.4;
stopa = Table[
{Opacity[v],Disk[xyG[T+dt*v],b]},{v,0.,1,1/24}];

stopa2 = Table[
{Opacity[v],Disk[xyG2[T+dt*v],b]},{v,0.,1,1/24}];
stopa3 = Table[
{Opacity[v],Disk[xyG3[T+dt*v],b]},{v,0.,1,1/24}];


Graphics[Join[
{Red}  ,stopa ,
{Blue} , stopa2,
{Darker@Green} , stopa3,
{Opacity[1],Yellow,Circle[{0,0},a]}
],
Background->Black,PlotRange->{{-a-b,a+b},{-a-b,a+b}}
],{T,0,60,0.1}];
ListAnimate[movie]

Export[NotebookDirectory[]<> "chaos.avi",movie]