# How to plot $(x,y)$ coordinate of projectile motion (with air-resistance)? [closed]

I am trying to plot a graph to show the difference in projectile motion when it has air-resistance and when it doesn't have air-resistance.

I set the mass, $v_0$ and $\theta$ as constant

I can plot the motion by using formulas.

$$x=v_0t\cos \theta,$$

$$y=v_0t\sin \theta -\frac12 g t^2.$$

How do I plot the projectile motion with air resistance?

• A simple Google search would suffice to find trajectory under air resistance. Look at Book or Wiki or Article Feb 16, 2014 at 7:02
• This seems to belong more in scicomp.stackexchange.com than it does here. Feb 16, 2014 at 10:31
• To emphasize what @PPG said, question on the programming aspects and the pure math of computational solutions are off topic. The physics part is obviously allowed, but in this case you have already done all that. Feb 16, 2014 at 15:44
• This question appears to be off-topic because it is about programming. Feb 16, 2014 at 15:44
• I will ask the mods on Computational Science if they want this (in my mind it is pretty elementary for that site). Feb 16, 2014 at 15:45

You use for example python and mathplolib like in the following (slightly modified) example taken from http://faculty1.coloradocollege.edu/~sburns/toolbox/ODE_II.html. You need a working python/pylab installation to run it via python trajectory.py (if you save the file with this name) from the command line.

The idea is to solve newtons equation of motion numerically using an ode-solver.

Choose for example a coordinate system with $y$-axis pointing upwards and $x$-axis to the left.

Then the $x$ and $y$-components of force on the particle with velocity $\vec{v}$ would be:

$F_x = -kv_x|\vec{v}|$ and $F_y = -g -kv_y |\vec{v}|$

where $k = \frac{1}{2}\rho_{\mathrm{air}} c_W A$ and $\rho_{\mathrm{air}}$ is the density of the air, $c_W$ is the drag coefficient and $A$ ist the area of the object perpendicular to the velocity vector.

The corresponding differential equation (from $\vec{F} = m \vec{a} = m \ddot{\vec{x}}$ is:

\begin{align} \ddot{x}(t)&= -kv_x(t)|\vec{v}(t)| \\ \ddot{y}(t) &= -g -kv_y(t) |\vec{v}(t)| \end{align}

Since most ode-solvers available can only handle first oder ode's you have to convert it to a first oder equation:

\begin{align*} v_x(t) &= \dot{x}(t)\\ v_y(t) &= \dot{y}(t)\\ \dot{v_x}(t)&= -kv_x(t)|\vec{v}(t)| \\ \dot{v_y}(t) &= -g -kv_y(t) |\vec{v}(t)| \end{align*}

"""
trajectory.py
"""
from pylab import *
from scipy.integrate import odeint

## set initial conditions and parameters
g = 9.81            # acceleration due to gravity
th = 45.            # set launch angle
th = th * pi/180.   # convert launch angle to radians
v0 = 10.0           # set speed
k = 0.3             # controls strength of air drag

x0=0                # specify initial conditions
y0=10
vx0 = v0*sin(th)
vy0 = v0*cos(th)

## define function to compute f(X,t)
def f_func(state,time):
f = zeros(4)    # create array to hold f vector
f[0] = state[2] # f[0] = x component of velocity
f[1] = state[3] # f[1] = x component of velocity
f[2] = - k*(f[0]**2 + f[1]**2)**(0.5)*f[0]         # f[2] = acceleration in x direction
f[3] = -g - k*(f[0]**2 + f[1]**2)**(0.5)*f[1]       # f[3] = acceleration in y direction
return f

## set initial state vector and time array
X0 = [ x0, y0, vx0, vy0]        # set initial state of the system
t0 = 0.
# tf = input("Enter final time: ")
# tau = input("Enter time step: ")
tf = 3
tau = 0.05

# create time array starting at t0, ending at tf with a spacing tau
t = arange(t0,tf,tau)

## solve ODE using odeint
X = odeint(f_func,X0,t) # returns an 2-dimensional array with the
# first index specifying the time and the
# second index specifying the component of
# the state vector

# putting ':' as an index specifies all of the elements for
# that index so x, y, vx, and vy are arrays at times specified
# in the time array
x = X[:,0]
y = X[:,1]
vx = X[:,2]
vy = X[:,3]

## plot the trajectory
figure(1)
clf()

# for the following two lines you need a working LaTeX environment, if you don't want LaTeX labels just comment the lines out
rc('text', usetex=True)
rc('font', family='serif')

plot(x,y)
xlabel('x')
ylabel('y')
savefig('output.png')
show()


Output (with LaTeX enabled):

For a more elementary approach, I would recommend the video analysis program tracker which is able to do the same and additionally overlay it to an experimental video of a projectile motion with air drag. To do so you need to create a new "dynamic particle model" and type in the components of force as explained above. The advantage is that you don't need to convert the second order differential equation to a first order equation since tracker does this for you internally already.

On the other hand to get publication ready quality I would have a look at some LaTeX based tools for example the pst-ode package or asymptote.

Basicly you can use any math-environment which allows ode-solving and plotting you want, for example commercial ones like mathematica, maple or matlab or free ones like scilab, octave, sage, maxima or the python libraries scipy and mathplotlib used above. The idea is always the same, you will have just a different syntax.

Edit:

Here is a different approach using python together with LaTeX and pgfplots (which is very powerful in customizing the axis styles etc.):

You need to run it with pdflatex --shell-escape plot_air_drag.tex if you save the file as plot_air_drag.tex:

\documentclass{article}

\usepackage{pgfplots}
\usepackage{python}

\pagestyle{empty}

\begin{document}

\begin{python}
from pylab import *
from scipy.integrate import odeint
from functools import partial

## set initial conditions and parameters
g = 9.81            # acceleration due to gravity
th = 45.            # set launch angle
th = th * np.pi/180.   # convert launch angle to radians
v0 = 20.0           # set speed

x0=0                # specify initial conditions
y0=0
vx0 = v0*sin(th)
vy0 = v0*cos(th)

## define function to compute f(X,t)
def f_func_k(k,state,time):
f = zeros(4)    # create array to hold f vector
f[0] = state[2] # f[0] = x component of velocity
f[1] = state[3] # f[1] = x component of velocity
f[2] = - k*(f[0]**2 + f[1]**2)**(0.5)*f[0]         # f[2] = acceleration in x direction
f[3] = -g - k*(f[0]**2 + f[1]**2)**(0.5)*f[1]       # f[3] = acceleration in y direction
return f

def gen_plot_data(s,k):
f_func = partial(f_func_k,k)

## set initial state vector and time array
X0 = [ x0, y0, vx0, vy0]        # set initial state of the system
t0 = 0.
# tf = input("Enter final time: ")
# tau = input("Enter time step: ")
tf = 3
tau = 0.05

# create time array starting at t0, ending at tf with a spacing tau
t = arange(t0,tf,tau)

## solve ODE using odeint
X = odeint(f_func,X0,t) # returns an 2-dimensional array with the
# first index specifying the time and the
# second index specifying the component of
# the state vector

# putting ':' as an index specifies all of the elements for
# that index so x, y, vx, and vy are arrays at times specified
# in the time array
x = X[:,0]
y = X[:,1]
vx = X[:,2]
vy = X[:,3]

xy = zip(x,y)
np.savetxt(s,xy,fmt='%0.5f')

gen_plot_data("plotdata1.dat",0.04)
gen_plot_data("plotdata2.dat",0)
\end{python}

\begin{tikzpicture}
\begin{axis}[no marks,samples=100,axis lines=center,xlabel=$x$,ylabel=$y$,enlargelimits]
\end{axis}
\end{tikzpicture}

\end{document}


Output:

To get a smoother plot you can change the tau variable in the python code and the samples variable in the latex-code to appropriate values.

• Thanks a lot. If I want to have graph look like this riotstories.co.uk/wp-content/themes/child2011/images/… how should i do sir Feb 16, 2014 at 9:30
• See my edit above Feb 16, 2014 at 10:16
• @Barbiyong To make the question clearer, you could the picture in the link into your original question! Feb 16, 2014 at 10:43
• If this problem has been given as homework (which I did this semester), the instructor almost certainly wants the students to perform the quadrature themselves (i.e. not to call a solver) because there are going to be following assignment on how and why the quadrature fails and how to correct them. Feb 16, 2014 at 15:49
• I am a bit confused about how to do the quadrature. It is clear that the vertical case (one dimension) can be separated and solved analytically. However with an arbitrary $\theta$ the problem becomes pretty complicated. You can solve it numerically as shown above, you can also use perturbation methods to get approximate solutions you can. The general problem is discussed in arxiv.org/pdf/1305.1283.pdf ... Feb 16, 2014 at 16:49

Try adding the deceleration caused due to air resistance in the $x$-direction equation. Remember that the air resistance is velocity dependent, so you'll have something like this:

$$F_{drag} = \frac{1}{2} \rho C_d A v^2$$

Which you would then put in the x-motion equation and integrate from time = zero to t