← all posts
MathematicsLinear AlgebraProofMachine Learning

Investigating Projectile Motion

Targeting solutions become significantly more difficult when drag is taken into account. We want to calculate a firing solution for any two position vectors with drag as a factor.

We want to find a way to calculate a solution without the usage of the As we will see, this requires the Lambert WW function.

Reference Frame

We consider two bodies: AA (the launcher) and BB (the target). Placing AA at the origin of our reference frame:

P⃗A=⟨0,0⟩,P⃗B=⟨xB, yB⟩.P_A = 0, 0 , P_B = x_B,\, y_B .

The projectile is launched from P⃗AP_A with initial velocity vector

v⃗0=v0cos⁡θ x^+v0sin⁡θ y^,v_0 = v_0\,x + v_0\,y,
where v0=∣v⃗0∣v_0 = |v_0| is the launch speed and θ is the elevation angle.

Solving an Ideal Projectile

For an ideal projectile, the differential equations below model velocity, where g=9.81 m/s2g = 9.81 m/s^2:

dvxdt=0dv_xdt = 0
dvydt=−g⋅dtdv_ydt = -g dt

The only rate of change in velocity is across the y-axis, thus the velocity along the x-axis remains constant.

Therefore, we can model vpx(t)v_px(t) as

vpx(t)=vpxiv_px(t) = v_pxi

To solve for vpy(t)v_py(t), we will solve the differential equation.

∫dvydt dt=∫−g dtvpy(t)=−gt+Cvpy(0)=vpyiC=vpyivpy(t)=−gt+vpyiarrayc dv_ydt \,dt = -g\,dt \\ v_py(t)= -gt + C \\ v_py(0)= v_pyi \\ C = v_pyi \\ v_py(t)= -gt + v_pyi array

Thus, we have solved our set of differential equations for velocity.

Now, given the elementary physics concept:

v(t)=drdtv(t)= drdt

We can relate our equations

vpy(t)=dpydtv_py(t)= dp_ydt
vpx(t)=dpxdtv_px(t)= dp_xdt

Creating a parametric equation to solve for this equation is ideal, as a value of x will directly correspond to a y value which we can check for displacement.

dpydt⋅dtdpx=dpydpxdp_ydt dtdp_x = dp_ydp_x
dpydpx=−g⋅dt+vpyivpxidp_ydp_x = -g dt + v_pyiv_pxi

However, dtdt is still in the equation and we need a relation between dtdt and pxp_x. Thus, we can integrate dpxdtdp_xdt

px(t)=∫dpxdt dt=∫vpxi dt=vpxi⋅t+Cp_x(t) = dp_xdt\, dt = v_pxi\, dt = v_pxi t + C

We can infer C=0C = 0 since the initial displacement of the projectile is 0. So we are left with

px(t)=vpxi⋅tp_x(t) = v_pxi t

Thus, from this equation, we can derive

dpxvpxi=dtdp_xv_pxi = dt

We substitute back into the differential equation

dpydpx=−g(dpxvpxi)+vpyivpxidp_ydp_x = -g(dp_xv_pxi) + v_pyiv_pxi

Which we simplify to

dpydpx=−g⋅dpxvpxi2+vpyivpxidp_ydp_x = -g dp_xv_pxi^2 + v_pyiv_pxi

We integrate the equation with respect to pxp_x

py(px)=∫dpydpx dpx=−g⋅px22vpxi2+vpyivpxi⋅px+Cp_y(p_x) = dp_ydp_x\, dp_x = -g p_x^22v_pxi^2 + v_pyiv_pxi p_x + C

Again, when px=0p_x = 0, py(px)=0p_y(p_x) = 0, thus C=0C = 0

py(px)=−g⋅px22vpxi2+vpyivpxi⋅pxp_y(p_x) = -g p_x^22v_pxi^2 + v_pyiv_pxi p_x

For a static target, we set py(px)=0p_y(p_x) = 0 and px=axp_x = a_x.

For simplicity, we will initially assume the displacement along the y-axis is 0.

0=axvpxi⋅(−g⋅ax2vpxi+vpyi)0 = a_xv_pxi (-g a_x2v_pxi + v_pyi)
0=−g⋅ax2vpxi+vpyi0 = -g a_x2v_pxi + v_pyi
−vpyi=−g⋅ax2vpxi-v_pyi = -g a_x2v_pxi
−2⋅vpyi⋅vpxi=−g⋅ax-2 v_pyi v_pxi = -g a_x
2⋅vpyi⋅vpxi=g⋅ax2 v_pyi v_pxi = g a_x

We need to solve for θ, thus we substitute the trigonometric identities:

vpxi=cos⁡(θ)⋅vpiv_pxi = () v_pi, vpyi=sin⁡(θ)⋅vpi v_pyi = () v_pi

2vpi2⋅sin⁡(θ)⋅cos⁡(θ)=g⋅ax2 v_pi^2 () () = g a_x
sin⁡(2θ)=g⋅axvpi2(2) = g a_xv_pi^2
θ=sin⁡−1 ⁣(g⋅axvpi2)2 = ^-1\!(g a_xv_pi^2)2

However, this will not consider displacement along the y-axis. Now, we will solve the more complex equation:

ay=g⋅ax22vpxi2+vpyivpxi⋅axa_y = g a_x^22v_pxi^2 + v_pyiv_pxi a_x
ay=g⋅ax22vpi2⋅cos⁡2(θ)+tan⁡(θ)⋅axa_y = g a_x^22v_pi^2 ^2() + () a_x
ay=g⋅ax22vpi2⋅(1+tan⁡2(θ))+tan⁡(θ)⋅axa_y = g a_x^22v_pi^2 (1 + ^2()) + () a_x
ay=g⋅ax22vpi2⋅tan⁡2(θ)+tan⁡(θ)⋅ax+g⋅ax22vpi2a_y = g a_x^22v_pi^2 ^2() + () a_x + g a_x^22v_pi^2
0=g⋅ax22vpi2⋅tan⁡2(θ)+tan⁡(θ)⋅ax+(g⋅ax22vpi2−ay)0 = g a_x^22v_pi^2 ^2() + () a_x + (g a_x^22v_pi^2 - a_y)

We use the quadratic formula to solve for the desired value of tan⁡(θ)()

tan⁡(θ)=−ax±ax2−4(g⋅ax22vpi2)(g⋅ax22vpi2−ay)g⋅ax2vpi2() = -a_x a_x^2 - 4(g a_x^22v_pi^2)(g a_x^22v_pi^2 - a_y)g a_x^2v_pi^2
tan⁡(θ)=vpi2±vpi4−g(g⋅ax2+2⋅vpi2⋅ay)g⋅ax() = v_pi^2 v_pi^4 - g(g a_x^2 + 2 v_pi^2 a_y)g a_x
θ=arctan⁡ ⁣(vpi2±vpi4−g(g⋅ax2+2⋅vpi2⋅ay)g⋅ax) = \!(v_pi^2 v_pi^4 - g(g a_x^2 + 2 v_pi^2 a_y)g a_x)

Thus, the following equation will solve for any stationary target given a displacement axa_x and aya_y. Solving for a moving target leads to a quartic equation, to be addressed in a future section.

Drag Definition

In this system, drag will be acting on this projectile where ζ=(0,1] = (0,1], where at ζ value of 1 is an ideal system with no drag. Drag is modeled by the differential:

dv⃗=v⃗⋅ζ⋅dtdv = v dt

We can rewrite this differential

dv⃗=dp⃗dt⋅ζ⋅dt=dp⃗⋅ζdv = dpdt dt = dp

We know that when the projectile's velocity is equal to 0, it has come to rest with respect to a certain axis. As the gravitational force is exclusive to the y-axis, we can model the maximum possible distance to be traveled along a certain axis to be.

v0x⃗ζ=Δpx⃗maxv_0_x = _x_max

We consider the two forces to be acting on the projectile to be the gravitational and the drag force, such that

a⃗x=−v⃗x⋅ζa⃗y=−g+−v⃗y⋅ζa_x = -v_x a_y = -g + -v_y

X-Axis Differential

Let us first derive a closed function for position along the x axis by taking the integral of the differential equation.

v⃗x(t)=∫dv⃗xdt dt=−v⃗x⋅ζ=e−ζt⋅v⃗x0v_x(t) = dv_xdt \ dt = -v_x = e^- t v_x_0

We integrate once more with respect to time to get a closed function for position with respect to the x-axis

p⃗x(t)=∫v⃗x(t) dt=−v⃗x0ζ⋅e−ζt+v⃗x0ζp_x(t) = v_x(t) \ dt = -v_x_0 e^- t + v_x_0

Which can be rewritten as:

p⃗x(t)=v⃗x0ζ⋅(1−e−ζt)p_x(t) = v_x_0 (1 - e^- t)

Differential with Gravity

We now need to integrate with respect to the y-axis, where gravity is a factor in this differential. Let g=10m/s2g = 10 m/s^2.

v⃗y(t)=∫a⃗ydt dt=(vy0⃗+gζ)⋅e−dx+g−vy0⃗ζv_y(t) = a_ydt \ dt = (v_y_0+g) e^-dx+g-v_y_0

We integrate once more to get a closed function for position with respect to the y-axis

p⃗y(t)=∫v⃗ydt dt=(v⃗y0+gζ)ζ(1 −e−ζt)+g−v⃗y0ζ⋅tp_y(t) = v_ydt \ dt = (v_y_0+g)(1\ -e^- t) + g-v_y_0 t

Finding a solution with respect to this axis is now incredibly difficult because we have an equation in the form of x=xm+qx = x^m + q

We will need to figure out a "trick" or a method to work around this.