Computing Assignment 2

Question 1

Here we need to find the linear path taken by \(A\) to intercept \(B\). We define angle \(\phi\) made by \(\boldsymbol{v}_A\) with the axis \(Ox\) as shown in the figure below and find \[ \boldsymbol{v}_A = v_A (-\cos\phi \boldsymbol{\hat{\imath}}+\sin\phi \boldsymbol{\hat{\jmath}})= \dot{x}_A \boldsymbol{\hat{\imath}}+ \dot{y}_A \boldsymbol{\hat{\jmath}} \] We integrate to find \[ x_A (t)= D - v_A t \cos\phi, \qquad y_A (t)= v_A t \sin\phi \] The position of \(B\) is given by \[ x_B (t)= 0, \qquad y_B (t) = v_B t \] \(A\) intercepts \(B\) at time \(t=T_1\) given by \[ v_A T_1 \cos\phi= D, \qquad v_A T_1 \sin\phi = v_B T_1 \] This gives 2 equations to find \((\phi,T_1)\): \[ \boxed{\sin\phi = \frac{v_B}{v_A}}, \qquad \boxed{T_1 = \frac{D}{\sqrt{v_A^2 -v_B^2}}} \] A solution exists only if \(\boxed{v_B < v_A}\). This is the necessary condition for \(A\) to catch up with \(B\).

Question 2.1. To find the ODEs which govern the motion of \(A\), we first find the unit vector which points from \(A (x_A,y_A)\) to \(B (0,v_B t)\): \[ \boldsymbol{\hat{e}}_{AB}= \frac{-x_A \boldsymbol{\hat{\imath}}+(v_B t- y_A) \boldsymbol{\hat{\jmath}}}{\sqrt{x_A^2 +(v_B t -y_A)^2}} \] Then we write the velocity of \(A\) as \(\boldsymbol{v}_A = \dot{x}_A \boldsymbol{\hat{\imath}}+\dot{y}_A \boldsymbol{\hat{\jmath}}= v_A \boldsymbol{\hat{e}}_{AB}\): this gives 2 equations \[ \dot{x}_A = - v_A \frac{x_A}{\sqrt{x_A^2 +(v_B t -y_A)^2}} = F (x_A,y_A, t) \qquad[1] \] and \[ \dot{y}_A = v_A \frac{(v_B t- y_A) }{\sqrt{x_A^2 +(v_B t -y_A)^2}} = G (x_A, y_A, t) \qquad [2] \]

2.2. If we denote \(x_n= x_A (t_n)\) and \(y_n = y_A (t_n)\) the values reached by \(x_A\) and \(y_A\) at time \(t=t_n= n \delta t\), Euler’s scheme gives the following iteration rules \[ x_{n+1} = x_{n} + \delta t F(x_n,y_n,t_n) \qquad (3) \] \[ y_{n+1} = y_{n} + \delta t G(x_n,y_n,t_n) \qquad (4) \] starting with the initial values \(x_0 = D\) and \(y_0=0\) at \(t_0=0\).

Given values of parameters \((D, v_A, v_B)\), we can create arrays of \(t\), \(x_A\) and \(y_A\) values by using the iteration equations (3-4) from \(t_0=0\) to \(t_N= N \delta t\) (parameter Tmax in the code). An estimation of parameter Tmax can be found by using time \(T_1\) (when \(v_A > v_B\)). To determine the time and corresponding position at which \(A\) catches up with \(B\), we need a stopping criterion: to accomplish this, we compute the distance \(|AB|\) at each time step, and stop the iterations as soon as \(|AB|< \epsilon\) is reached. \(\epsilon\ll 1\) is a prescribed small parameter (parameter epsilon in the code). The smaller \(\epsilon\) is chosen, the smaller \(\delta t\) needs to be.

To determine the array of position values of \(A\) according to equations [1-2], we create a function which returns \(F\) and \(G\) given values of \(x,y\) and \(t\).

def Velocity(x,y,t):
    # returns an velocity components F and G as functions of position (x,y) and time t
    X = -x
    Y = vB*t-y
    dist= np.sqrt(X^2 +Y^2)
    F= vA*X/dist
    G= vA*Y/dist
    return F, G

To determine the acceleration of \(A\) (and its maximum value) along the path, we use the fact that \(A\) moves at constant speed, and hence \(|\boldsymbol{a}_A| = v_A^2 /\rho\). To estimate the radius of curvature \(\rho\) at a point \((x_n,y_n)\) of the trajectory, we compute the radius of the circle which passes through \((x_n,y_n)\) and its nearest neighbors \((x_{n-1},y_{n-1})\) and \((x_{n+1},y_{n+1})\). If \(\delta t\) is very small, this circle is very good approximation of the osculating circle.

def FindCircle(x1, y1, x2, y2, x3, y3) :
# returns the radius of the circle passing through  
# three given points
    x12 = x1 - x2  
    x13 = x1 - x3  
    y12 = y1 - y2  
    y13 = y1 - y3
    y31 = y3 - y1  
    y21 = y2 - y1  
    x31 = x3 - x1  
    x21 = x2 - x1  
    sx13 = x1**2 - x3**2  
    sy13 = y1**2 - y3**2  
    sx21 = x2**2 - x1**2  
    sy21 = y2**2 - y1**2  
    f = (sx13*x12+sy13*x12+sx21*x13+sy21*x13)/(2*(y31*x12-y21*x13))
    g = (sx13*y12+sy13*y12+sx21*y13+sy21*y13)/(2*(x31*y12-x21*y13))  
    c = -x1**2 - y1**2 - 2 * g * x1 - 2 * f * y1
  
    # eqn of circle is x^2 + y^2 + 2*g*x + 2*f*y + c = 0  
    # where center is (h = -g, k = -f) and  
    # radius r as r^2 = h^2 + k^2 - c  
    h = -g;  
    k = -f;  
    r = np.sqrt(h * h + k * k - c)  
    return r

The full script is available here.

2.3 Results.