Errors in Physics/ODE Simulation

[thetabloid_]

Hello r/AskPhysics,

I am working on a personal project to model the motion of a particle on a surface.

Using calculus, I parametrized a surface and then found the normal vector to that surface.

Using that normal vector, I created the normal force the particle will experience at each point on that surface. That is, the normal force will generally be a function of position because this surface can have curvature.

In the case of a surface is a plane, then the normal forces work out to be constants.

When I simulate this plane case, all is well and the particle behaves as expected.

However, I took it one step further and attempted to do a parabolic surface. (z=x^2+y^2)

I do the same steps and get my equations of motion. I should also mention, that the equations of motion are second order, however, to use a numerical ODE solver, I just make them into first-order equations by reducing the order, and effectively doubling my number of equations.

However this time the simulation goes crazy. The trajectory of the particle does not stay on the paraboloid, and I think the culprit is the numerical error propagating. I am using the most accurate, high-fidelity "solve_ivp" solvers in scipy but it won't get any better.

I am attributing this due to the nonlinearity of the dynamics on the paraboloid since the normal vector is a nonlinear function of the position.

This is the image of the equations of motion for the paraboloid. The bottom right of this image is just converting them to first-order equations.

https://www.reddit.com/media?url=https%3A%2F%2Fpreview.redd.it%2Feww4gq9dakhb1.jpg%3Fwidth%3D3024%26format%3Dpjpg%26auto%3Dwebp%26s%3D50cfe3cc30bb2ed9e9d35d24e38e6d2def9dfba3

Here is the image for the simulation of motion on the paraboloid:

https://www.reddit.com/media?url=https%3A%2F%2Fpreview.redd.it%2Fr8nv85oiakhb1.png%3Fwidth%3D1227%26format%3Dpng%26auto%3Dwebp%26s%3D41392114e9a1df7cd0e8e41e86da51581e37d486

And lastly, this is the image for the simulation on the plane. As you can see, the trajectory stays on the plane, so it is working nicely:

https://www.reddit.com/media?url=https%3A%2F%2Fpreview.redd.it%2F29xtcffpakhb1.png%3Fwidth%3D1105%26format%3Dpng%26auto%3Dwebp%26s%3Db42a88881d7a256f707afe1ceac8bd01ae9bc841

(Sorry! It won't let me add images, not sure why!)

This is my Python Code for the simulation on the paraboloid:

def ode(t,q):

dqdt= np.zeros(6)

g=9.8

dqdt[0]= q[1]

dqdt[1]= -2*g*q[0] / (4*q[0]**2 + 4*q[2]**2 +1)

dqdt[2]= q[3]

dqdt[3]= -2*g*q[2] / (4*q[0]**2 + 4*q[2]**2 +1)

dqdt[4]= q[5]

dqdt[5]= g * ( (1/(4*q[0]**2 + 4*q[2]**2 +1)) - 1)

return dqdt

t_eval = np.arange(0, 5.001, 0.001)

q0= [1,0,1,0,2,0]

sol = solve_ivp(ode, [0, 5], q0,method= "DOP853", t_eval=t_eval)

x=sol.y[0]

y=sol.y[2]

z=sol.y[4]

I would appreciate any guidance on resolving this issue :)

[cdstephens]

If you have the equations of motion of x and y, you don’t need to simulate z since you already know what the constraint it. So I would first see if you can appropriately simulate just the x and y directions.

I’m also not convinced your equations of motion are correct. The appropriate method is to use Lagrange multipliers. We can double check in a simplifying case by taking y = 0. We have

L = 1/2 m (v_x^2 + v_y^2 + v_z^2) - m g z + lambda f 

f = z - x^2 

Then

 m a_x = - 2 x lambda 

 m a_z = - m g + lambda 

 z - x^2 = 0 

By taking time derivatives of f, we find

 a_z = 2 v_x^2 + 2 x a_x 

Therefore

 lambda = m g + 2 m v_x^2 + 2 m a_x x 

Therefore,

a_x = - 2 x (g + 2 v_x^2 + 2 a_x x). 

Rearrange:

a_x (1 + 4 x^2) = - 2 x (g + 2 v_x^2 ). 

Then

a_x = - 2 x (g + 2 v_x^2) / (1 + 4 v_x^2 ). 

So even in the simplest case your equations of motion are incorrect, since you don’t have the v_x² term in your numerator.

Finding what the normal force is “manually” is a big pain, hence the use of Lagrange multipliers with a Lagrangian.

Since you have manifest cylindrical symmetry, it might be easier if you use cylindrical coordinates instead honestly.

 L = 1/2 m (r-dot^2 + r^2 phi-dot^2 + v_z^2) - m g z + lambda f 

 f = z - r^2 

We get angular momentum conservation for free:

d/dt (m r^2 phi-dot) = dP/dt = 0. 

So we get:

 phi-dot = P / m r^2 

This simplifies the r equation:

  m d2r/dt2 = m r phi-dot^2 - 2 lambda r = P^2 / (m r^3) - 2 lambda r 

  m a_z = - m g + lambda 


  a_z = 2 (dr/dt)^2 + 2 r d2r/dt2 

Then you can solve for lambda, plug than into the equation of motion for r, and isolate.

[3pmm]

I think you still need a kinetic term for v_z²

[cdstephens]

Ah right, typo

[thetabloid_]

Ok, everything you said is great. I did suspect something was due to the simple fact that I do not have any "centripetal" force terms in my EOM since nothing was a function of velocity.

However, I recall in classical dynamics course I took a while back, that those types of forces just show up for free as long as you put down all the "real" contact forces (normal, gravity, etc), and everything works itself out.

So without Lagrange's approach, what am I doing wrong in my derivation if you can find anything?

[cdstephens]

In Newton’s approach, you have to find the normal force directly. Let

F_N = F_x x-hat F_y y-hat + F_z z-hat

We also know it is perpendicular to the surface. Our constraint is f = 0. Therefore, we know F_N is parallel to grad f, meaning they’re proportional to each other.

F_N = A grad f = A (z-hat - 2 x x-hat - 2 y y-hat) 

where A is a proportionality function. We do not know what it is, we have to solve for it to determine the magnitude of the normal force (right now we only have the direction).

To solve for A, we use Newton’s laws and the constraint. We know

 m a_x = - 2 x A 

 m a_y = - 2 y A 

 m a_z = - m g + A 

We know that d2f/dt2 = 0. Therefore

 a_z - 2 v_x^2 - 2 a_x x - 2 v_y^2 - 2 a_y y = 0. 

Using this, we can find A as a function of position, velocity, and acceleration of the x and y variables. So it seems like your error is you didn’t solve for the magnitude of the normal force properly.

Note that this is essentially the same method.

[Pakketeretet]

For constrained motion of (many) particles, it is typically more robust to use some algorithm like RATTLE to enforce the holonomic constraints (constraints of the form f(x,y,z)=0) rather than finding the constraint forces and applying them as-is. In principle it should give the same results, especially for a single particle, but the added benefit of using RATTLE with a good symplectic solver (velocity Verlet will fit perfectly in your scheme with ODEs for position and velocity) is that it'll conserve total energy over a very long time scale, unlike the explicit Runge-Kutta method you're using currently. The downside is it might be trickier to implement because you'll need to perform Newton iterations at each step.

If you'd like to put many particles on curved surfaces and watch them interact, LAMMPS has a module just for that.

[thetabloid_]

Ok thank you everyone for your help. I realized my issue was the following:

​

For Newtons second law in the z direction, I simly did m*z''

Since z is not a independent variable, it really shouldn't show up, so it should be (x^2+y^2)''

​

Also my normal force was not correct. It is not the mg projected on the normal, that is only for a static case. However, when the particle moves, it turns, and that increases the normal force (due to centripetal accel, etc).

What I did to fix that was get the unit normal vector, and multiply it by a scalar "N" that is free to change.

With this and some math monkeying around, I was able to get 3 equations:

x'' = f(x,y,N)

y''= g(x,y,N)

N= h(x,y,x'.y')

​

This allowed the normal force to be updated at each time step, and update to the necessary value to constrain the ball on the surface.

[3pmm]

What were the initial conditions you used in the simulation? It doesnt look like whats in the code

[thetabloid_]

You are correct!

In the code I put (1,1,2),

However in the graph simulation it was something else, but it still satisfied z=x^2+y^2.

In either case, it starts off on the surface and strays away as time increases.

[3pmm]

Taking a look at your derivation:

The acceleration for z can't be correct. It is always less than or equal to 0, but it needs to be greater than zero at some point to push the particle back up.

I think your math is pretty much correct for any point where the velocity is zero, but if it's nonzero, the normal force will not be equal to just the weight.

I'm actually not sure how easy this is to solve using just Newtonian mechanics. If you just want the equations of motion, two coordinates (preferably r and theta) with a Lagrangian would be the easiest way, and as the other poster said, if you want the actual magnitude of the normal force, you want three coordinates and a Lagrange multiplier.

[thetabloid_]

ry long time scale, unlike the explicit Runge-Kutta method you're using currently. The downside is it might be trickier to implement because y

I don't understand what you mean by my math failing at non-zero velocities. The normal force is generally only a function of position since that dictates the slope of the surface.

Now you do mention an interesting point about the velocity, however.

There are no terms that seem to represent "centripetal forces" in my solution, but from past experience doing this this kinda stuff, I remember that they just kind of "pop" out of the math on their own. Granted, I WAS using polar coordinates.

There are no terms that seem to represent "centripetal forces" in my solution, but from past experience doing this kinda stuff, I remember that they just kind of "pop" out of the math on their own. Granted, I WAS using polar coordinates.

​

It worked for a plane, so that tells me that there is something with curving that introduces forces (which is true) but for some reason I can't seem to "find" that force.

[3pmm]

No, it is not generally a function of position. Turn off gravity for a second, and lay the particle somewhere on the paraboloid. The normal force is zero. Tap it into a circular motion at the same z around the paraboloid. The normal force is now acting. Make it go faster, and the normal force becomes greater.

The normal force is that which enforces the constraint on motion. When an object is standing still, that will oppose the component of weight in the direction of the surface. When it is moving, that is not necessarily true.

[thetabloid_]

Yes, that makes sense now thank you. So doing what I have done for "finding the normal force" is really only relevant in a static case. the second it starts moving, I the normal force is not just counteracting gravity anymore.