# Plot3D Equation of motion of a particle in a surface

I'm new in this forum. This semester I've been trying to learn how to use Mathematica for my Classical Mechanics course, but I'm still a novice concerning this software. In particular today I was doing a problem about a water drop on a surface $$z=x^2-y^2$$ and I came with Lagrangian Mechanics to the following equations:

$$(1+4x^2)\ddot{x}-4y\ddot{y}x+4x\dot{x}-4\dot{y}^2x-2gx=0$$ $$(1+4y^2)\ddot{y}-4x\ddot{x}y+4y\dot{y}-4\dot{x}^2y+2gx=0$$

where $$g=9,8$$ and arbitrary initial conditions. I tried a lot of codes, but I always had mistakes. Please, can you help me how to Plot3D this trajectory? Thanks in advance.

• g = 9.8; sol = NDSolve[ {(1 + 4 x[t]^2) x''[t] - 4 y[t] y''[t] x[t] + 4 x[t] x'[t] - 4 y'[t]^2 x[t] - 2 g x[t] == 0, (1 + 4 y[t]^2) y''[t] - 4 x[t] x''[t] y[t] + 4 y[t] y'[t] - 4 x'[t]^2 y[t] + 2 g x[t] == 0, x[0] == 1, x'[0] == 1, y[0] == 1, y'[0] == 1} , {x, y}, {t, 0, 5}]; Plot[x[t] /. First@sol, {t, 0, 5}]; Plot[y[t] /. First@sol, {t, 0, 5}]; ParametricPlot[{x[t], y[t]} /. First@sol, {t, 0, 5}] ? Commented Oct 30, 2020 at 9:04
• @cvgmt ah we seem to have done the same thing! Commented Oct 30, 2020 at 9:08
• @chris I don't know how to give the initial conditions :) Commented Oct 30, 2020 at 9:15
• People here generally like users to post code as Mathematica code instead of just images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. You may find this meta Q&A helpful Commented Oct 30, 2020 at 12:30

If I understand correctly, you have a point mass on a surface $$z=x^2-y^2$$ with gravitation. But then your equations are wrong.

For simplicity we choose $$m=1$$, then

$$L= 1/2 (x'^2+y'^2)-g (x^2-y^2)$$

And this gives the ODE:

$$x''+2g x=0$$ and $$y'' - 2 g y =0$$

In MMA:

g = 9.81; sol =
NDSolve[{x''[t] + 2 g x[t] == 0, y''[t] - 2 g y[t] == 0, x[0] == 1,
x'[0] == 1, y[0] == 0, y'[0] == 0.1}, {x, y}, {t, 0, 1}];

Show[ParametricPlot3D[{x[t], y[t], x[t]^2 - y[t]^2} /. First@sol, {t,
0, 1}, AxesLabel -> {"x", "y", "z"}]
, ParametricPlot3D[{x, y, x^2 - y^2}, {x, -1, 4}, {y, -1, 4}]]


Or with different initial conditions:

g = 9.81; sol =
NDSolve[{x''[t] + 2 g x[t] == 0, y''[t] - 2 g y[t] == 0, x[0] == 1,
x'[0] == 1, y[0] == -1, y'[0] == 4.5}, {x, y}, {t, 0, 1}];

Show[ParametricPlot3D[{x[t], y[t], x[t]^2 - y[t]^2} /. First@sol, {t,
0, 1}, AxesLabel -> {"x", "y", "z"}]
, ParametricPlot3D[{x, y, x^2 - y^2}, {x, -1, 4}, {y, -1, 4}]]


Take care with initial conditions, because the surface is pretty steep in y direction and the mass point speeds easily away.

• You reversed engineered the question? :) Nice answer anyway! Commented Oct 30, 2020 at 16:30

Something like this ?

eqn = {(1 + 4 x[t]^2) D[x[t], t, t] - 4 y[t] D[y[t], t, t] x[t] +
4 x[t] D[x[t], t] -
4 x[t] D[y[t], t]^2 - 2 g x[t] == 0,
(1 + 4 y[t]^2) D[y[t], t, t] - 4 x[t] D[x[t], t, t] y[t] +
4 y[t] D[y[t], t] -
4 y[t] D[x[t], t]^2 + 2 g x[t] == 0
} /. g -> 10


eqn2 = Solve[eqn, {x''[t], y''[t]}] // FullSimplify //
First // # /. Rule -> Equal &;


Then

A[t_] = NDSolveValue[
Flatten[{eqn2,
x[0] == 0, x'[0] == 1,
y[0] == 1/10, y'[0] == 1/10}],
{x[t], y[t]}, {t, 0, 0.5}];
ParametricPlot[A[t], {t, 0, 0.5}]


Note that the equation is stiff so it starts to do something wrong eventually.

You can explore different initial conditions as follows

A[t_] = Table[NDSolveValue[
Flatten[{eqn2,
x[0] == 0, x'[0] == 1,
y[0] == i/20, y'[0] == 1/10}],
{x[t], y[t]}, {t, 0, 0.5}], {i, 1, 5}];

ParametricPlot[A[t] // Evaluate, {t, 0, 0.5},
PlotStyle -> NestList[Lighter, Red, 10]]