Zadanie polega na wyznaczeniu toru kulki (punktu materialnego) poruszajacej
się bez tarcia po zadanej powierzchni z=f[x,y] w jednorodnym polu grawitacyjnym g.

Dowolna powierzchnia

z  =  f[x[t], y[t]]

Paraboloida obrotowa

z = -(1 + x[t]^2 + y[t]^2)

Pewna bardzo pofaldowana powierzchnia

z = FresnelS[x[t]^2] FresnelC[y[t]^2]

Inna przykładowa powierzchnia

z = Exp[-x[t]^2 - y[t]^2] Cos[x[t]]

Płaszczyzna prostopadła do siły grawitacyjnej

z  =  0

Dowolna płaszczyzna nachylona pod pewnym kątem

z =  A x[t]  +  B y[t]  +  C

UWAGA: wyżej należy wybrać tylko JEDNĄ z podanych funkcji z[x[t],y[t]]

Tu sprawdzamy czy MATHEMATICA poprawnie stosuje regułę łańcuchową:

D[z, t] - D[z, x[t]] x'[t] - D[z, y[t]] y'[t]

Simplify[%]

Układ równań opisujących tor znajdziemy z równań Lagrange-Eulera

Potrzebujemy wyrażenia na energię kinetyczną T kulki równą T=1/2 m v^2 wyrażoną przez położenie kulki x[t] i y[t]

T =  1/2 m (x'[t]^2 + y'[t]^2 + D[z, t]^2)

Wyrażenie na energię potencjalną V jest proste:

V  =  -g z

Funkcja Lagrange'a to różnica energii kinetycznej i potencjalnej wyrażona przez x[t], y[t], x'[t] i y'[t]

L = T - V

Teraz tworzymy układ równań Lagrange-Eulera.

Lewa strona pierwszego r. L-E

eq1L  = D[D[L, x'[t]], t]

Prawa strona pierwszego r.L-E

eq1R =  D[L, x[t]]

Lewa strona drugiego r.L-E

eq2L  = D[D[L, y'[t]], t]

Prawa strona drugiego r.L-E

eq2R =  D[L, y[t]]

Pierwsze z równań L-E. Podstawiamy masę kulki m=1 i przyspieszenie ziemskie
g=1.  Można oczywiście uzyć bardziej realistycznych wartości ;-)

eq11  =  (eq1L == eq1R)  /.  m -> 1  /.  
        g -> 1

Drugie z równań:

eq22  =  (eq2L == eq2R)  /.  m -> 1  /.  
        g -> 1

Układ 2 równań Lagrange-Eulera który musimy rozwiazać aby znaleźć tor kulki:

sys  =  {eq11, eq22}

Warunki początkowe.Podajemy położenie i prędkość w chwili t=0

ic = {x[0] == 1.5,y[0] == 1.0,x'[0]==-0.01,y'[0]==0.5}

W nielicznych przypadkach (np. ruch po płaszczyźnie) możliwe jest znalezienie rozw. analitycznego:

DSolve[sys,  {x[t], y[t]}, t]

W wiekszości przypadków równania różniczkowe zwyczajne można rozwiązać wyłącznie numerycznie. MATHEMATICA robi to bardzo dobrze automatycznie dobierając algorytm i kontrolując dokładność. Na ogół nie potrzeba używać żadnych dodatkowych opcji. Składnia NDSolve jest bardzo podobna do DSolve. Podstawowe różnice to: 1) warunki początkowe ( ic)  muszą być podane 2) wszystkie parametry (w tym war. początkowe!) muszą mieć nadane konkretne wartości numeryczne (w przykładzie masa i g) 3) musimy podać konkretny przedział czasowy  ( t)w którym ma zostać znalezione rozwiazanie (0.0, 4.0).

tr  =  NDSolve[Union[sys, ic], {x[t], y[t]}, {t, 0.0, 15}]

Tu możemy zobaczyć tor w rzucie z góry:

ParametricPlot[{x[t], y[t]}  /.  tr, {t, 0, 15}]

A tutaj w 3D:

h =  (z  /.  tr)[[1]]

Liczba 0.01 ma na celu narysowanie później toru nad powierzchnią (inaczej jej nie widać)

d  =  ({x[t], y[t], -h + 0.01}  /.  tr)[[1]]

tr1  =  
      ParametricPlot3D[d,  {t, 0, 15}, PlotPoints -> 1000, 
        ViewPoint -> {1.0, 2.0, 2.0}]

Oto powierzchnia po której porusza się kulka:

pow  =  
      Plot3D[-z, {x[t], 1.2, 2.2}, {y[t], 0.7, 1.2}, Mesh -> False, 
        PlotPoints -> 50, ViewPoint -> {1.0, 2.0, 2.0}]

A tu kulka na powierzchni:

Show[pow, tr1]

Zachęcam do eksperymentowania z różnymi wartościami początkowymi  i powierzchniami