Computing Assignment 1

Question 1

We found in Lecture 1 that the trajectory connecting \(A(x_A,0)\) to \(B(0,h)\) with launch speed/angle \((v_0,\alpha)\) must satisfy a quadratic equation in \(X=\tan\alpha\): \[ C X^2 +B X +A = 0, \qquad A= 2h + \frac{gx_A^2}{v_0^2}, \quad B= 2x_A , \quad C= \frac{gx_A^2}{v_0^2} \]
This quadratic equation has in general two roots unless \(B^2 -4AC <0\). This last condition shows that a path can be found only if \(v_0\) exceeds a minimum value: \[ v_0 > v_{0min} = \sqrt{gh} (1+ (1+x_A^2/h^2)^{1/2}) \] Given values of parameters \(x_A\) and \(h\), no path can be found for \(v_0 < v_{0min}\).

Python code: For given values of \((x_A,v_0,h)\) we compute the values of angle \(\alpha\) by solving the quadratic equation in \(X=\tan\alpha\) which guarantees that \(A\) and B are connected. Since 3 values of \(v_0\) are given, we set up an array for \(v_0\). We add a fourth value which corresponds to \(v_{0min}\), below which there is no solution. The fifth value of \(v_0\) is the optimal path. On a graph, we will plot the 5 trajectories and the rectangular obstacle. It is convenient to set up a function which computes the path for given values of \((x_A,v_0,\alpha)\): this functions returns 2 arrays of \((x,y)\) values associated with this path.

def Projectile(X, v, alpha):
    # returns an array of [x, y] values from t=0 to tfinal
    # starting at A(X,0) with initial speed v and launch angle
    # alpha
    t0 = 0
    tfinal = -X/(v*np.cos(alpha))  # final time to reach point B(0,h) starting from
                                    # point A(xA,0)
    t = np.linspace(t0, tfinal, Npt)    
    x =  X + v*t*np.cos(alpha)
    y =  v*t*np.sin(alpha)-0.5*g*t**2
    return [x, y]

For each value of \(v_0\) we solve the quadratic equation for angle \(\alpha\), call function Projectile and use the values returned by this function to plot the path.

We obtain this plot (the yellow path is the optimal path found in Q.2):

Question 2: Q.1 shows that, for given values of \((x_A,h)\), there is a minimum launch speed v0min necessary to reach point \(B\). Also we note that v0min decreases as a function of xA. There is good reason to believe that there will be a minimum value of xA below which the paths requires higher initial speed v0.

Based on these observations, we can build a scheme with two nested loops:

First we create a loop of xA values, we increase xA as the loop proceeds (that is, the position of \(A\) is closer to \(O\)).
Inside this loop, we create an array of v0 values from v0min to v0max.
We then create a v0 loop, by decreasing v0 as this second loop proceeds: Inside this loop, we determine the paths connecting \(A(x_A,0)\) to \(B(0,h)\) by finding the two launch angles. These paths are either acceptable or not.
- A path is acceptable if no point \(P(x,y)\) of the path is inside the region defined by the building: \(D< x <0 , \qquad 0< y <h\). For the rectangular building, it is straightforward to test: we find the time t_hit when the value \(x=-d\) is reached. Then we find the quantity y_hit as the value of \(y\) at \(t=\) t_hit. If y_hit > h, the path is acceptable. There is no need to call function Projectile.
- If the path is acceptable, we test if the corresponding v0 value is lower that our previously recorded minimum. If it is, this is our new minimum.

Eventually, as the computations proceed along these 2 loops, there are no longer smaller values of v0. We can get closer to the obstacle, but it takes a larger speed v0 to reach \(B\). The scheme converges very fast to the solution. You can refine your solution by adding more values to the loops.

This gives the following script:

for xA in X:
    v0min=  np.sqrt(g*h*(1+np.sqrt(1.0+xA**2/h**2)))
    v0max= 1.2*v0min
    Vvalues= np.linspace(-v0max, -v0min, Nvvalues, endpoint=False)# array of v0 values
    Vvalues = -Vvalues

    for V0 in reversed(Vvalues):    # v0 loop from v0max to v0min, excluding v0min
                                    # (at v0min, discrim is a very small number, which can
                                    # be <0 due to machine precision)

        A =2*h + g*xA**2/V0**2
        B = 2*xA
        C= g*xA**2/V0**2
        discrim=B**2-4*A*C

        alpha1 = np.arctan(0.5*(-B+np.sqrt(discrim))/C)  # launch angle 1 
        t1_hit= -(d+xA)/np.cos(alpha1)/V0         # corresponding time to reach x=-d
        y1_hit= -0.5*g*t1_hit**2 + V0*t1_hit*np.sin(alpha1)  # height at x=-d

        alpha2 = np.arctan(0.5*(-B-np.sqrt(discrim))/C) # launch angle 2
        t2_hit= -(d+xA)/np.cos(alpha2)/V0    # time to reach x=-d
        y2_hit= -0.5*g*t2_hit**2 + V0*t2_hit*np.sin(alpha2) # height at x=-d
        
        if (y1_hit > h): # test if path #1 is acceptable
            if(V0 < Vmin):  # if acceptable, test if V0 is smaller than Vmin
                Vmin=V0     # if smaller, then this is the new minimum
                Xmin=xA      # record the corresponding value of xA
                AlphaMin= alpha1 # record the corresponding value of alpha
                print('Rectangle Case: xA=',xA,'[m], v0=',V0,'[m/s], alpha=',alpha1*180/np.pi,'degrees')

        if (y2_hit > h):  # repeat for the second path
            if(Vmin>V0):
                Vmin=V0
                Xmin=xA
                AlphaMin=alpha2
                print('Rectangle Case: xA=',xA,'[m], v0=',V0,'[m/s], alpha=',alpha2*180/np.pi,'degrees')

Numerically we find the optimal path is given by \[ \boxed{ x_A = -23.95 \text{ [m]}, \quad v_0= 27.29 \text{ [m/s]}, \quad \alpha= 75.04 \text{ degrees} } \]

The second case proceeds the same way, except for the testing of the paths. To test if a path is acceptable, we need to guarantee that no point \((x,y)\) of the path joining \(A\) to \(B\) satisfies the condition \[ x^2 + (y-h/2)^2 < (h/2)^2 \] In this case, for each value of \((x_A, v_0)\) in the nested loops, we need to make a call to function Projectile and test that the returned \((x,y)\) values are acceptable or not. To do this, we use a function check:

def check(DistanceArray, value): 
    return(all(x > value for x in DistanceArray))

where the array DistanceArray contains the values of \(x^2+ (y-h/2)^2\) and the call to check is made with value = (h/2)**2. If the returned value is True, the path is acceptable. This case is more costly to run than the first case. We find the following values for the optimal path \[ \boxed{ x_A = -24.30 \text{ [m]}, \quad v_0= 26.98 \text{ [m/s]}, \quad \alpha= 72.94 \text{ degrees} } \] and the optimal path given by:

The full script is available here.