Why is my simulation of projectile trajectory with air resistance wrong?

[fisicist]

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:

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()

Result:

This shouldn't happen... As you can see, the curves bend upward, i.e. are not concave. Even without knowing the exact solution $y(x)$, it is easy to prove it should be a concave function.
The model is $\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$
By the chain rule, $\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)$
Therefore, at least as long as $\dot{y} \geq 0$, 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: $\frac{d^2 y}{d^2 x} (x(t)) = -\frac{1}{2 \dot{x}^2} \leq 0$

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?

 

[Dale]

Therefore, at least as long as $\dot{y} \geq 0$, this should be a concave function.

Hmm, I am not seeing that, where do you get that condition from?

 

[fisicist]

Hmm, I am not seeing that, where do you get that condition from?

$1/\dot{x}^2 \geq 0$ manifestly, so it remains to check that $-2 c \dot{y} \sqrt{\dot{x}^2+\dot{y}^2}-\frac 1 2 \leq 0$ in order for $d^2 y/ d x^2 \leq 0$ (i.e. concavity). But $-2 c \sqrt{\dot{x}^2 +\dot{y}^2}$ is always negative, so if $\dot{y} \geq 0$, everything is negative.

 

[Dale]

What about the 1/2?

 

[A.T.]

Hmm, I am not seeing that, where do you get that condition from?

An informal argument for concavity: The drag is always tangential so it doesn't change the direction, while gravity is always bending the path down.

Looks like some lift sneaked in somewhere.

 

[fisicist]

Hey sorry, I've made a sign error. The two friction terms compensate each other, not add up. Actually not surprising: Without gravity, motion is along a straight line.
I will correct it. Correct is
$d^2 y / dx^2 = - \frac{1}{2 \dot{x}^2}$
Still, it's concave. In fact, it is *always* concave, even if I throw the stone downwards.

 

[A.T.]

As you can see, the curves bend upward, i.e. are not concave.

I would call bending upward and then downward "concave". Bending only downward would be "convex", and what you would expect here.

Have you tried coding a simple Euler integration yourself?

 

[fisicist]

I would call bending upward and then downward "concave". Bending only downward would be "convex", and what you would expect here.

Have you tried coding a simple Euler integration yourself?

The standard definition used in mathematics (in fields such as calculus of variations, 'convex analysis' or 'convex optimization') is that a function is convex if its epigraph $\operatorname{epi}f := \lbrace (x, y) : y \geq f(x) \rbrace$ is a convex set. For a twice differentiable function, this is equivalent to $f'' \geq 0$. A function $f$ is concave if $-f$ is convex.

 

[A.T.]

The standard definition used in mathematics (in fields such as calculus of variations, 'convex analysis' or 'convex optimization') is that a function is convex if its epigraph $\operatorname{epi}f := \lbrace (x, y) : y \geq f(x) \rbrace$ is a convex set. For a twice differentiable function, this is equivalent to $f'' \geq 0$. A function $f$ is concave if $-f$ is convex.

Ok, thanks. I would suggest coding a simple Euler integration. It should never bend upwards.

 

[fisicist]

Thank you to everyone who replied here!
I found the error. It's not a mathematical or numerical curiosity, but a simple coding error. Namely,

np.absolute(vel)

does not give $\vert \dot{\mathbf x} = \sqrt{\dot{x}^2+\dot{y}^2}$, as I thought it would, but $(\vert \dot x \vert, \vert \dot y \vert )^T$, and

np.absolute(vel) * vel

does not give $\dot{\mathbf{x}} \vert \dot{\mathbf{x}} \vert$, as I thought it would, but $(\dot{x} \vert \dot{x} \vert, \dot{y} \vert \dot{y} \vert )^T$.

The correct way to do it is

np.sqrt(vel.dot(vel)) * vel

 

[William Crawford]

Thank you to everyone who replied here!
I found the error. It's not a mathematical or numerical curiosity, but a simple coding error. Namely,

np.absolute(vel)

does not give $\vert \dot{\mathbf x} = \sqrt{\dot{x}^2+\dot{y}^2}$, as I thought it would, but $(\vert \dot x \vert, \vert \dot y \vert )^T$, and

np.absolute(vel) * vel

does not give $\dot{\mathbf{x}} \vert \dot{\mathbf{x}} \vert$, as I thought it would, but $(\dot{x} \vert \dot{x} \vert, \dot{y} \vert \dot{y} \vert )^T$.

The correct way to do it is

np.sqrt(vel.dot(vel)) * vel

Congratulation! So what is the ideal angle when you include air resistance?

 

[fisicist]

For anyone interested in the physics: Qualitatively, the result is similar.

Qualitatively, the best angle decreases rapidly the higher the friction is:

 

[Mister T]

Run your sim on a baseball and see if 40 deg gives the maximum range.

 

[A.T.]

The correct way to do it is

np.sqrt(vel.dot(vel)) * vel

Or:

np.linalg.norm(vel) * vel

https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.norm.html

For anyone interested in the physics: Qualitatively, the result is similar.
View attachment 247485
Qualitatively, the best angle decreases rapidly the higher the friction is:
View attachment 247487

Regarding your stone throwing contest:

You could plot the optimal angle as function of the object's density for a given realistic object shape/size and throw speed. It could be that for stones the angle is pretty close to 45° at the speed you can throw, while for light balls (like tennis) it makes sense to aim lower.

Bio-mechanics complicates matters because you won't achieve the same throw speeds for any angle.

Finally, your friend might simply throw much faster, and thus further, even at sub optimal angles.