Exercise_06 : Enhanced version exercise 2.10

To hit the target with the minimum firing velocity

冉峰 弘毅班 2014301020064

---

Abstract

• With the method of Euler we can solve many problem using computer. A further discussion on the last question, we take the effect of wind and air density into consideration to calculate the minimum firing velocity to hit a target with different altitude and distance.

---

Background

From the Computational Physics we know the method of Euler.Now I want to take the Euler method into the projectile movement calculation:
Newton’s second law in two spatial dimensions :

$$\frac{d^2x}{dt^2}=\frac{F_x}{m}=\frac{F_{drag,x}}{m}\qquad\frac{d^2y}{dt^2}=\frac{F_y}{m}=-g+\frac{F_{drag,y}}{m}$$
Write each of these second-order equations as two firest-order differential equations:
$$\frac{dx}{dt}=v_x\qquad\frac{dv_x}{dt}=\frac{F_{drag,x}}{m}$$
$$\frac{dy}{dt}=v_y\qquad\frac{dv_y}{dt}=-g+\frac{F_{drag,y}}{m}$$
Finite difference form:
$$x_{i+1}=x_i+v_{x,i}\Delta t\qquad v_{x,i+1}=v_{x,i}+\frac{F_{drag,x}}{m}\Delta t$$
$$y_{i+1}=y_i+v_{y,i}\Delta t\qquad v_{y,i+1}=v_{y,i}-g\Delta t+\frac{F_{drag,y}}{m}\Delta t$$
For the $F_{drag}^*$ with the effect of the wind we know : $$F_{drag,x}^*=-B_2|\vec{v}-\vec{v}_{wind}|(v_x-v_{wind})$$
$$F_{drag,y}^*=-B_2|\vec{v}-\vec{v}_{wind}|v_y$$
Also we consider the air density with the adiabatic approximation: $$F_{drag}=\frac{\rho}{\rho_0}F_{drag}^*(y=0)=(1-\frac{ay}{T_0})^\alpha F_{drag}^*(y=0)$$
So we have:
$$x_{i+1}=x_i+v_{x,i}\Delta t\qquad v_{x,i+1}=v_{x,i}-\frac{B_2|\vec{v_i}-\vec{v}_{wind}|(v_{x,i}-v_{wind})}{m}(1-\frac{ay_i}{T_0})^\alpha\Delta t$$
$$y_{i+1}=y_i+v_{y,i}\Delta t\qquad v_{y,i+1}=v_{y,i}-g\Delta t+\frac{B_2|\vec{v_i}-\vec{v}_{wind}|v_{y,i}}{m}(1-\frac{ay_i}{T_0})^\alpha\Delta t$$
where $\alpha=2.5 \quad,\quad a=6.5\times 10^{-3}K/m$

---

The Main Body

• code 1

import math
xx = 1000 #x position of the target,unit:m
yy = 100 #y position of the target,unit:m
g = 9.8 
a = (6.5)*(math.pow(10,-3)) 
Bm = 4*math.pow(10,-5) #B2/m,unit:m^(-1)
T0 = 300 
vw = -20 #velocity of wind, unit:m/s
va = 1000 #initial_velocity,unit:m/s
sv = [va]
sa = [0]
s = xx
for i in range(21):
    v = va 
    v0 = 0
    angle = 25 + i*0.001 #initial_angle=25°,add 2° each time
    for j in range(20):
        x = [0] #initial_x,unit:m
        y = [0] # initial_y,unit:m
        dt = 0.1 #time_step,unit:s
        vx = [v*math.cos(math.pi*angle/180)]
        vy = [v*math.sin(math.pi*angle/180)]
        while (x[-1]<=xx): 
            if y[-1] >100000000:
                c = 0
            else:
                c = math.pow((1-(a*y[-1]/T0)),2.5)
            Fx = -Bm*math.sqrt((vx[-1]-vw)*(vx[-1]-vw)+vy[-1]*vy[-1])*(vx[-1]-vw)*c
            Fy = -Bm*math.sqrt((vx[-1]-vw)*(vx[-1]-vw)+vy[-1]*vy[-1])*vy[-1]*c
            vx.append(vx[-1]+Fx*dt)
            vy.append(vy[-1]+(-g+Fy)*dt)  
            x.append(x[-1]+vx[-2]*dt)
            y.append(y[-1]+vy[-2]*dt)
        r = (xx-x[-2])/(x[-1]-xx)
        y[-1] = (y[-2]+r*y[-1])/(r+1)
        x[-1] = xx
        if y[-1]>yy:
            vt = v
            v = v0+(vt-v0)/2
        if y[-1]<yy:
            v0 = v
            v = v0+(vt-v0)/2
    sv.append(v)
    sa.append(angle)
vmin = min(sv)
l = sv.index(vmin)
amin = sa[l]
print('The x position of the target',xx)
print('The y position of the target',yy)
print('velocity of wind',vw)
print ('The minimum initial velocity:'+str(vmin)+'m/s')
print('angle:'+str(amin)+'°')

Result

• result

The x position of the target 1000
The y position of the target 100
velocity of wind -20
The minimum initial velocity:129.1036605834961m/s
angle:25.02°

---

The x position of the target 1000
The y position of the target -100
velocity of wind -20
The minimum initial velocity:104.19178009033203m/s
angle:25.019°

---

The x position of the target 100
The y position of the target -100
velocity of wind -20
The minimum initial velocity:20.12920379638672m/s
angle:25.0°

---

The x position of the target 100
The y position of the target 100
velocity of wind -20
The minimum initial velocity:20.131091947405366m/s
angle:25.0°

---

• code 2

import math
import pylab as pl
import math
import time
class cannon:
    def __init__(self, vel_0=1000, time_of_duration=1, vel_0_step=0.01, time_step=0.01, m=100, B=0.01, vel_wind=2, theta=0, step=0.01, value_of_x=0, value_of_y=0, g=10, T0=300):
        self.v0 = [vel_0]
        self.theta = [theta]
        self.vx = [self.v0[0]]
        self.vy = [self.v0[0]]
        self.x = [value_of_x]
        self.y = [value_of_y]
        self.t = [0]
        self.g = g
        self.dt = time_step
        self.time = time_of_duration
        self.nsteps = int(time_of_duration // time_step + 1)
        self.dtheta = step
        self.nsteps2 = int(90//step + 1)
        self.B = B
        self.m = m
        self.dv = vel_0_step
        self.vw = vel_wind
        self.T0 = T0
        self.n = 0
    def calculate(self):
        for k in range(100):
            tmpv0 = self.v0[k] + self.dv
            self.v0.append(tmpv0)
            self.theta = [0]
            self.vx = [0]
            self.vy = [0]
            for j in range(self.nsteps2):
                tmptheta = self.theta[j] + self.dtheta
                self.theta.append(tmptheta)
                self.x = [0]
                self.y = [0]
                self.vx[0] = self.v0[k] * math.cos(self.theta[j] * math.pi / 180)
                self.vy[0] = self.v0[k] * math.sin(self.theta[j] * math.pi / 180)
                for i in range(self.nsteps):
                    tmpvx = self.vx[i] - self.dt * pow((1-0.0065 * self.y[i]/self.T0),2.5) * self.B * math.sqrt(float(self.vx[i]-self.vw) * float(self.vx[i]-self.vw) + self.vy[i] * self.vy[i]) * (self.vx[i] - self.vw)/self.m
                    tmpvy = self.vy[i] - self.g * self.dt + self.dt * pow((1-0.0065 * self.y[i]/self.T0),2.5) * self.B * \
                    math.sqrt(float(self.vx[i]-self.vw) * float(self.vx[i]-self.vw)+ self.vy[i] * self.vy[i]) * self.vy[i]/self.m
                    tmpx=self.x[i] + self.vx[i] * self.dt
                    tmpy=self.y[i] + self.vy[i] * self.dt
                    self.vx.append(tmpvx)
                    self.vy.append(tmpvy)
                    self.x.append(tmpx)
                    self.y.append(tmpy)
                    self.t.append(self.t[i] + self.dt)
                    if int(self.x[i]) == 10:
                        if self.y[i] > 10:
                            self.n = 1
                            break
                if self.n == 1:
                    break
            if self.n == 1:
                self.show_results()
                print("The minimum of velocity is:%f", self.v0[k])
                break
    def show_results(self):
        pl.plot(self.x, self.y)
        pl.xlabel('x ($m$)')
        pl.ylabel('y ($m$)')
        pl.show()
a = cannon()
a.calculate()

result

The x position of the target 10000
The y position of the target 10000
velocity of wind 2
The minimum initial velocity:1000.10999m/s

---

But the code have a problem which is too slowly .

Conclusion

• With the help the of method of Euler, we can solve many physics problem by computer easy.In my opinion, we can solve almost all the problem we meet in calculating the cannon trajectories. What the most important is that we can skilled use the python to solve the problem. Maybe it's the significance of we learning computationalphysics.

---

Thanks For

• Lesson plan of Cai Hao:Chapter 2
• xiongpeiyu-2014301200262