The trajectory of cannon shell

air drag density trajectory

Abstract

The Euler method we used to treat the bicycle problem can easily be generalized to deal with motion in two spatial dimensions.To be specific,we consider a projectile such as a shell shot by a cannon.

Introduction

Generally, when we ignore the effects of air drag, the trajectory of a firing cannon ball is a parabola obviously. However, if we consider air drag and even the diverse density of air at different altitudes, the situation will become much more complicated than the previous one.

Text

If we ignore air resisttance,the equations of motion can be obtained from Newton's second law,can be written as
$\begin{cases}\frac{d^2x}{dt^2}=0\\\frac{d^2y}{dt^2}=-g\end{cases}$
we can use Euler methods to solve this second-order ODE by writing them as followings
$\begin{cases}\frac{dx}{dt}=v_x\\\frac{dv_x}{dt}=o\\\frac{dy}{dt}=v_y\\\frac{dv_y}{dt}=-g\end{cases}$

import math  
import pylab as pl
class bicycle:
    def __init__(self, power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/6),\
       initial_Vy=50*math.sin(math.pi/6)):
        self.v = [initial_velocity]
        self.x=[initial_x]
        self.y=[initial_y]
        self.Vx=[initial_Vx]
        self.Vy=[initial_Vy]
        self.t = [0]
        self.m = mass
        self.p = power
        self.dt = time_step
    def run(self):
        while(self.y[-1]>=0):
            self.x.append(self.x[-1]+self.Vx[0]*self.dt)
            self.y.append(self.y[-1]+self.Vy[-1]*self.dt)
            self.Vy.append(self.Vy[-1]-9.8*self.dt)
        print("yn=",  self.y[-1])
        print("y(n-1)=",self.y[-2])
        print("xn=",self.x[-1])
        print("x(n-1)=",self.x[-2])
    def show_results(self):
        pl.plot(self.x, self.y)
        pl.title('projectile motion without air resistance')
        pl.xlabel('x')
        pl.ylabel('y')
        pl.show()
a = bicycle()
a.run()
a.show_results()

let $r=-y_{n-1}/y_{n}$
$x_l=\frac{x_{n-1}+rx_{n}}{r+1}$

If change the initial angles of projectile motion,we can get the following results.

import numpy as np
import math  
import pylab as pl
class bicycle:
    def __init__(self, power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/6),\
       initial_Vy=50*math.sin(math.pi/6)):
        self.v = [initial_velocity]
        self.x=[initial_x]
        self.y=[initial_y]
        self.Vx=[initial_Vx]
        self.Vy=[initial_Vy]
        self.t = [0]
        self.m = mass
        self.p = power
        self.dt = time_step
    def run(self):
        while(self.y[-1]>=0):
            self.x.append(self.x[-1]+self.Vx[0]*self.dt)
            self.y.append(self.y[-1]+self.Vy[-1]*self.dt)
            self.Vy.append(self.Vy[-1]-9.8*self.dt)
        print("yn=",  self.y[-1])
        print("y(n-1)=",self.y[-2])
        print("xn=",self.x[-1])
        print("x(n-1)=",self.x[-2])
class diff_step_check(bicycle):    
    def show_result1(self,style='b'):
        pl.plot(self.x, self.y,style,label='30')
        pl.title('projectile motion with air resistance')
        pl.xlabel('x')
        pl.ylabel('y')
        pl.legend(loc='best')
from pylab import *
a=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/6),\
       initial_Vy=50*math.sin(math.pi/6))
a.run()
a.show_result1('b')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='35')
        pl.title('projectile motion with air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/6*(7/6)),\
       initial_Vy=50*math.sin(math.pi/6*(7/6)))
b.run()
b.show_result2('g')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='40')
        pl.title('projectile motion with air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
c=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/6*(4/3)),\
       initial_Vy=50*math.sin(math.pi/6*(4/3)))
c.run()
c.show_result2('r')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='45')
        pl.title('projectile motion with air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
d=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/4),\
       initial_Vy=50*math.sin(math.pi/4))
d.run()
d.show_result2('c')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='50')
        pl.title('projectile motion with air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/4*(10/9)),\
       initial_Vy=50*math.sin(math.pi/4*(10/9)))
b.run()
b.show_result2('m')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='55')
        pl.title('projectile motion with air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/4*(11/9)),\
       initial_Vy=50*math.sin(math.pi/4*(11/9)))
b.run()
b.show_result2('y')
class diff_step_check(bicycle):
    def show_result3(self,style='b'):
        pl.plot(self.x, self.y,style,label='60')
        pl.title('projectile motion with air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
c=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/3),\
       initial_Vy=50*math.sin(math.pi/3))
c.run()
c.show_result3('k')
show()

Then we add the air drag $F_{drag}=-B_2v^2$,
the results is as followings

import numpy as np
import math  
import pylab as pl
class bicycle:
    def __init__(self, power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/6),\
       initial_Vy=50*math.sin(math.pi/6)):
        self.v = [initial_velocity]
        self.x=[initial_x]
        self.y=[initial_y]
        self.Vx=[initial_Vx]
        self.Vy=[initial_Vy]
        self.t = [0]
        self.m = mass
        self.p = power
        self.dt = time_step
    def run(self):
        while(self.y[-1]>=0):
            self.x.append(self.x[-1]+self.Vx[-1]*self.dt)
            self.Vx.append(self.Vx[-1]-self.v[-1]*self.Vx[-1]/self.m*self.dt)
            self.y.append(self.y[-1]+self.Vy[-1]*self.dt)
            self.Vy.append(self.Vy[-1]-9.8*self.dt-self.v[-1]*self.Vy[-1]/self.m*self.dt)
        print("yn=",  self.y[-1])
        print("y(n-1)=",self.y[-2])
        print("xn=",self.x[-1])
        print("x(n-1)=",self.x[-2])
class diff_step_check(bicycle):    
    def show_result1(self,style='b'):
        pl.plot(self.x, self.y,style,label='30')
        pl.title('projectile motion with air resistance')
        pl.xlabel('x')
        pl.ylabel('y')
        pl.legend(loc='best')
from pylab import *
a=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/6),\
       initial_Vy=50*math.sin(math.pi/6))
a.run()
a.show_result1('b')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='35')
        pl.title('projectile motion with air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/6*(7/6)),\
       initial_Vy=50*math.sin(math.pi/6*(7/6)))
b.run()
b.show_result2('g')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='40')
        pl.title('projectile motion with air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
c=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/6*(4/3)),\
       initial_Vy=50*math.sin(math.pi/6*(4/3)))
c.run()
c.show_result2('r')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='45')
        pl.title('projectile motion with air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
d=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/4),\
       initial_Vy=50*math.sin(math.pi/4))
d.run()
d.show_result2('c')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='50')
        pl.title('projectile motion with air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/4*(10/9)),\
       initial_Vy=50*math.sin(math.pi/4*(10/9)))
b.run()
b.show_result2('m')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='55')
        pl.title('projectile motion with air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/4*(11/9)),\
       initial_Vy=50*math.sin(math.pi/4*(11/9)))
b.run()
b.show_result2('y')
class diff_step_check(bicycle):
    def show_result3(self,style='b'):
        pl.plot(self.x, self.y,style,label='60')
        pl.title('projectile motion with air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
c=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/3),\
       initial_Vy=50*math.sin(math.pi/3))
c.run()
c.show_result3('k')
show()

problem 2.9-- the cannon shell trajectory

Firstly,we add the air drag.

import numpy as np
import math  
import pylab as pl
class bicycle:
    def __init__(self, power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/6),\
       initial_Vy=700*math.sin(math.pi/6)):
        self.v = [initial_velocity]
        self.x=[initial_x]
        self.y=[initial_y]
        self.Vx=[initial_Vx]
        self.Vy=[initial_Vy]
        self.t = [0]
        self.m = mass
        self.p = power
        self.dt = time_step
    def run(self):
        while(self.y[-1]>=0):
            self.x.append(self.x[-1]+self.Vx[-1]*self.dt)
            self.Vx.append(self.Vx[-1]-math.exp(-self.y[-1]/(10**4))*4/(10**5)*self.v[-1]*self.Vx[-1]/self.m*self.dt)
            self.y.append(self.y[-1]+self.Vy[-1]*self.dt)
            self.Vy.append(self.Vy[-1]-9.8*self.dt-math.exp(-self.y[-1]/(10**4))*4/(10**5)*self.v[-1]*self.Vy[-1]/self.m*self.dt)
            self.v.append(math.sqrt(self.Vx[-1]**2+self.Vx[-1]**2))
        print("yn=",  self.y[-1])
        print("y(n-1)=",self.y[-2])
        print("xn=",self.x[-1])
        print("x(n-1)=",self.x[-2])
class diff_step_check(bicycle):    
    def show_result1(self,style='b'):
        pl.plot(self.x, self.y,style,label='30')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x')
        pl.ylabel('y')
        pl.legend(loc='best')
from pylab import *
a=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/6),\
       initial_Vy=700*math.sin(math.pi/6))
a.run()
a.show_result1('b')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='35')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/6*(7/6)),\
       initial_Vy=700*math.sin(math.pi/6*(7/6)))
b.run()
b.show_result2('g')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='40')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
c=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/6*(4/3)),\
       initial_Vy=700*math.sin(math.pi/6*(4/3)))
c.run()
c.show_result2('r')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='45')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
d=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/4),\
       initial_Vy=700*math.sin(math.pi/4))
d.run()
d.show_result2('c')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='50')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/4*(10/9)),\
       initial_Vy=700*math.sin(math.pi/4*(10/9)))
b.run()
b.show_result2('m')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='55')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/4*(11/9)),\
       initial_Vy=700*math.sin(math.pi/4*(11/9)))
b.run()
b.show_result2('y')
class diff_step_check(bicycle):
    def show_result3(self,style='b'):
        pl.plot(self.x, self.y,style,label='60')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
c=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/3),\
       initial_Vy=700*math.sin(math.pi/3))
c.run()
c.show_result3('k')
show()

then we combine the air drag with the reduced air density$\rho=\rho_0exp(-y/y_0)$under the isothermal atmosphere

import numpy as np
import math  
import pylab as pl
class bicycle:
    def __init__(self, power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/6),\
       initial_Vy=700*math.sin(math.pi/6)):
        self.v = [initial_velocity]
        self.x=[initial_x]
        self.y=[initial_y]
        self.Vx=[initial_Vx]
        self.Vy=[initial_Vy]
        self.t = [0]
        self.m = mass
        self.p = power
        self.dt = time_step
    def run(self):
        while(self.y[-1]>=0):
            self.x.append(self.x[-1]+self.Vx[-1]*self.dt)
            self.Vx.append(self.Vx[-1]-math.exp(-self.y[-1]/(10**4))*4/(10**5)*self.v[-1]*self.Vx[-1]/self.m*self.dt)
            self.y.append(self.y[-1]+self.Vy[-1]*self.dt)
            self.Vy.append(self.Vy[-1]-9.8*self.dt-math.exp(-self.y[-1]/(10**4))*4/(10**5)*self.v[-1]*self.Vy[-1]/self.m*self.dt)
            self.v.append(math.sqrt(self.Vx[-1]**2+self.Vx[-1]**2))
        print("yn=",  self.y[-1])
        print("y(n-1)=",self.y[-2])
        print("xn=",self.x[-1])
        print("x(n-1)=",self.x[-2])
        r=-self.y[-2]/self.y[-1]
        xl=(self.x[-2]+r*self.x[-1])/(r+1)
        print("零点为",xl)
class diff_step_check(bicycle):    
    def show_result1(self,style='b'):
        pl.plot(self.x, self.y,style,label='30')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x')
        pl.ylabel('y')
        pl.legend(loc='best')
from pylab import *
a=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/6),\
       initial_Vy=700*math.sin(math.pi/6))
a.run()
a.show_result1('b')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='35')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/6*(7/6)),\
       initial_Vy=700*math.sin(math.pi/6*(7/6)))
b.run()
b.show_result2('g')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='40')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
c=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/6*(4/3)),\
       initial_Vy=700*math.sin(math.pi/6*(4/3)))
c.run()
c.show_result2('r')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='45')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
d=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/4),\
       initial_Vy=700*math.sin(math.pi/4))
d.run()
d.show_result2('c')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='50')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/4*(10/9)),\
       initial_Vy=700*math.sin(math.pi/4*(10/9)))
b.run()
b.show_result2('m')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='55')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/4*(11/9)),\
       initial_Vy=700*math.sin(math.pi/4*(11/9)))
b.run()
b.show_result2('y')
class diff_step_check(bicycle):
    def show_result3(self,style='b'):
        pl.plot(self.x, self.y,style,label='60')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
c=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/3),\
       initial_Vy=700*math.sin(math.pi/3))
c.run()
c.show_result3('k')
show()

we can also get the distance form approximate calculation.

It's easy to see that from $30^。$to $55^。$,the result is increasing,while from$55^。$to $60^。$,is downing.So the max angle is around $50^。$to $60^。$

Then change the angle by $2^。$ each time from $50^。$to $60^。$

the max angle is around $52^。$to $56^。$

Then change the angle by $1^。$ from $52^。$to $56^。$

so we can guass around the $55^。$ give the maximum range.

Thanks

By myself.