- Summary
- Simulation of trajectories with friction gives wrong results

Hey!

This started very harmless... A friend and I were throwing stones in a lake. Mine didn't get very far, he was teasing me "What was the ideal angle again?". Of course, I know it should be 45°. I replied in jest: "That's because I'm considering air resistance!" Then we had a discussion what the ideal angle should be with air resistance considered.

So I wrote a simulation in Python. It was just for fun... But I encountered an issue that I don't quite understand: The simulated trajectory is wrong. I hope someone can help me understand where this is coming from, and perhaps give me some hint how to fix it.

My code is:

Result:

This shouldn't happen... As you can see, the curves bend upward, i.e. are not concave. Even without knowing the exact solution [itex]y(x)[/itex], it is easy to prove it should be a concave function.

The model is [itex]\ddot{x} = - c \dot{x}\sqrt{\dot{x}^2+\dot{y}^2}, \ddot{y} = - c \dot{y} \sqrt{\dot{x}^2+\dot{y}^2} - \frac 1 2 [/itex]

~~By the chain rule, [itex]\frac{d^2 y}{d^2 x} (x(t)) = \frac{\ddot{y}}{\dot{x}^2} - \dot{y}\frac{\ddot{x}}{\dot{x}^3} = \frac{1}{\dot{x}^2} \left( -2 c \dot y \sqrt{\dot{x}^2+\dot{y}^2} - \frac 1 2 \right)[/itex]~~

Therefore, at least as long as [itex]\dot{y} \geq 0[/itex], this should be a concave function. If I through it downwards with high speed, it might get convex there, but not as long as the motion is upward.

CORRECTION: By the chain rule: [itex]\frac{d^2 y}{d^2 x} (x(t)) = -\frac{1}{2 \dot{x}^2} \leq 0[/itex]

So, what is wrong with the simulation? I know the time subdivision is not very high. But that shouldn't matter, because also in the discretized model, the velocity vector should be shortened and then get a downward contribution in every step. That's not what I see in the plot.

Any ideas, anyone?

This started very harmless... A friend and I were throwing stones in a lake. Mine didn't get very far, he was teasing me "What was the ideal angle again?". Of course, I know it should be 45°. I replied in jest: "That's because I'm considering air resistance!" Then we had a discussion what the ideal angle should be with air resistance considered.

So I wrote a simulation in Python. It was just for fun... But I encountered an issue that I don't quite understand: The simulated trajectory is wrong. I hope someone can help me understand where this is coming from, and perhaps give me some hint how to fix it.

My code is:

Python:

```
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
c = 7
v = 5
# f gives phase space velocity
# z = [x, y, v_x, v_y]
# f(t, z) = [v_x, v_y, a_x, a_y]
def f(t, z):
global c
vel = z[2:]
friction = -c * np.absolute(vel) * vel
gravity = [0, -0.5]
accel = friction + gravity
return np.concatenate((vel, accel))
def hit_ground(t, z): return z[1]
hit_ground.terminal = True
hit_ground.direction = -1
t_eval = np.linspace(0, 100, 10000)
def v_init(theta):
global v
return v * np.array([np.cos(theta), np.sin(theta)])
thetas = np.linspace(0.0, np.pi/2, 181)
dists = np.array([])
plt.subplot('211')
for counter,theta in enumerate(thetas):
init = np.concatenate((np.array([0., 0.]), v_init(theta)))
sol = solve_ivp(f, [0, 100], init, t_eval = t_eval, events=hit_ground)
if counter % 40 == 0:
plt.plot(sol.y[0,:], sol.y[1,:], label='{:4.0f} °'.format(theta*180/np.pi))
plt.gca().quiver(sol.y[0,::70], sol.y[1,::70], 1000*sol.y[2,::70], 1000*sol.y[3,::70])
dists = np.append(dists, [sol.y[0,-1]])
#plt.gca().set_aspect('equal', adjustable='box')
plt.legend(loc='upper left')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.subplot('212')
plt.plot(thetas*180/np.pi, dists, 'r-')
best_angle = thetas[np.argmax(dists)]
plt.gca().axvline(x = best_angle*180/ np.pi)
print('{:3.1f}'.format(best_angle*180/ np.pi))
plt.xlabel('angle [°]')
plt.ylabel('distance')
plt.show()
```

This shouldn't happen... As you can see, the curves bend upward, i.e. are not concave. Even without knowing the exact solution [itex]y(x)[/itex], it is easy to prove it should be a concave function.

The model is [itex]\ddot{x} = - c \dot{x}\sqrt{\dot{x}^2+\dot{y}^2}, \ddot{y} = - c \dot{y} \sqrt{\dot{x}^2+\dot{y}^2} - \frac 1 2 [/itex]

Therefore, at least as long as [itex]\dot{y} \geq 0[/itex], this should be a concave function. If I through it downwards with high speed, it might get convex there, but not as long as the motion is upward.

CORRECTION: By the chain rule: [itex]\frac{d^2 y}{d^2 x} (x(t)) = -\frac{1}{2 \dot{x}^2} \leq 0[/itex]

So, what is wrong with the simulation? I know the time subdivision is not very high. But that shouldn't matter, because also in the discretized model, the velocity vector should be shortened and then get a downward contribution in every step. That's not what I see in the plot.

Any ideas, anyone?

Last edited: