I have an object with an initial velocity in 2D space represented by a vector. I want to calculate a constant acceleration with a given magnitude required for the object to (potentially) change directions and pass through a given point (doesn't matter what angle it intersects the point or its speed when it does so, just that it intersects the point). I'm having a lot of trouble trying to figure out an equation or procedure to calculate this acceleration.
I'm running this in a physics engine, so the object has a position and a velocity (x speed and y speed). I am trying to find the direction of the acceleration required for the object to intersect the point. There are no other forces on the object. The velocity represents the distance the object travels in a single "step" of the engine and the acceleration represents the amount that the velocity changes in each step.
I tried to find a constant velocity by calculating the number of steps it would take for the object to reach the point on a single axis using the maximum acceleration, then normalizing the number of steps to get an acceleration in two dimensions. This is what I came up with:
total distance = distance covered by the initial velocity + distance gained by acceleration
so
$$p=vx+a*\frac{x(x+1)}{2}$$
where p = position, v = initial velocity, and x = step
I used this equation ($1+2+3+4+...=\frac{x(x+1)}{2}$) to calculate the distance gained by acceleration, since each new element can easily represent a step of the engine. I then used Wolfram Alpha to solve for x and got this:
$$x=-\frac{-\sqrt{a^2+8ap+4av+4v^2}+a+2v}{2a}$$
where $a\neq0$
This equation works both when the object's initial velocity is opposite from the target point and when its initial velocity is towards the point.
I had a lot of issues getting the object to precisely intersect the point and ultimately decided that this approach wasn't going to produce an accurate result. After that I came up with the idea that I could create some kind of parabolic curve based on the object's initial velocity and the target's position that would help me calculate the acceleration needed. Here's what I mean. I have no idea where I would even start for that, though.
Finally, I decided that it might be possible for me to use trigonometry to calculate the angles needed, but I don't really know where to start with that either. I don't know the distances that the object will have to travel on a single axis to intersect the point, since the distance changes based on the distance of the other axis (ie if the object travels faster on the x axis, it will travel slower on the y axis, and vice versa).
What is the best way to approach this problem?
A quick background refresher:
In continuous time, the velocity $v(t)$ and position $x(t)$ with a constant acceleration $a$ are given by $$\left\lbrace\begin{aligned} v(t) &= \int_{0}^{t} a \, d \tau = v_0 + a t \\ x(t) &= \int_{0}^{t} v(\tau) \, d \tau = x_0 + v_0 t + \frac{1}{2} a t^2 \end{aligned}\right.\tag{1}\label{NA1}$$
In discrete time, the velocity and position are given by difference equations. There are three most common definitions:
A. Position updated before velocity: $$\left\lbrace\begin{aligned} x_{i+1} &= x_{i} + v_{i} \\ v_{i+1} &= v_{i} + a \\ \end{aligned}\right. \implies x_N = x_0 + N v_0 + \frac{N (N - 1)}{2} a \tag{2a}\label{NA2a}$$
B. Velocity updated before position: $$\left\lbrace\begin{aligned} v_{i+1} &= v_{i} + a \\ x_{i+1} &= x_{i} + v_{i+1} \\ \end{aligned}\right. \implies x_N = x_0 + N v_0 + \frac{N(N+1)}{2} a \tag{2b}\label{NA2b}$$
C. Velocity updated before and after position: $$\left\lbrace\begin{aligned} v_{i+\frac{1}{2}} &= v_{i} + \frac{1}{2} a \\ x_{i+1} &= x_{i} + v_{i+\frac{1}{2}} \\ v_{i+1} &= v_{i+\frac{1}{2}} + \frac{1}{2} a \\ \end{aligned}\right. \implies x_N = x_0 + N v_0 + \frac{N^2}{2} a \tag{2c}\label{NA2c}$$
In this last case $\eqref{NA2c}$, velocity updates are essentially half a time step out of sync with the position updates, but the position function is most similar to the continuous time one, $\eqref{NA1}$.
(While the formulae above are shown in one dimension, you can use them in 2D or 3D just as well. The time steps ($i$ and $N$) are the same (shared) across all dimensions, but each dimension has their own acceleration $a$. The actual magnitude of the acceleration is then $\sqrt{a_x^2 + a_y^2}$ for 2D, and $\sqrt{a_x^2 + a_y^2 + a_z^2}$ for 3D.)
The problem at hand is to find the smallest $0 \le N \in \mathbb{N}$ and $a_x$ and $a_y$, when $x_0$, $y_0$, $v_{x0}$, $v_{y0}$ are known, but acceleration is limited, $a_x^2 + a_y^2 \le A_{max}^2$.
If we solve $\eqref{NA2a}$ for $a$ for both dimensions, we get $$\left\lbrace\begin{aligned} a_x &= \frac{2 (x_N - x_0 - N v_{x0})}{N (N - 1)} \\ a_y &= \frac{2 (y_N - y_0 - N v_{y0})}{N (N - 1)} \\ \end{aligned}\right . \tag{3a}\label{NA3a}$$ Solving $\eqref{NA2b}$ similarly yields $$\left\lbrace\begin{aligned} a_x &= \frac{2 (x_N - x_0 - N v_{x0})}{N (N + 1)} \\ a_y &= \frac{2 (y_N - y_0 - N v_{y0})}{N (N + 1)} \\ \end{aligned}\right . \tag{3b}\label{NA3b}$$ and $\eqref{NA2c}$ yields $$\left\lbrace\begin{aligned} a_x &= \frac{2 (x_N - x_0 - N v_{x0})}{N^2} \\ a_y &= \frac{2 (y_N - y_0 - N v_{y0})}{N^2} \\ \end{aligned}\right . \tag{3c}\label{NA3c}$$ To find out the approximate number of steps (which we need to round upwards to not exceed the specified acceleration $A_{max}$), we could try and solve $$a_x^2 + a_y^2 \le A_{max}^2 \tag{4}\label{NA4}$$ for $N \in \mathbb{R}$ (so the actual number of time steps needed would be rounded up, $\left\lceil N \right\rceil$).
Unfortunately, while there is an algebraic solution (it is a form of a fourth-degree function that does have algebraic solutions), it is too complicated to be useful.
Fortunately, we can use the bisection method (or binary search) to very efficiently find the smallest $N$, and therefore also the largest acceleration $(a_x , a_y)$ not exceeding $A_{max}$ in magnitude, that reaches the target coordinates in minimum number of time steps.
Note that in practice, square root and division operations take much more time than addition or multiplication. Because the bisection method/binary search will have to evaluate the same expression potentially dozens of times, it is a good idea to "optimize" $\eqref{NA4}$ first.
If we use $\Delta_x = x_N - x_0$ and $\Delta_y = y_N - y_0$, then $\eqref{NA4}$ optimizes to: $$(\Delta_x - N v_{x0})^2 + (\Delta_y - N v_{y0})^2 - N^2 (N-1)^2 \frac{A_{max}^2}{4} \le 0 \tag{5a}\label{NA5a}$$ $$(\Delta_x - N v_{x0})^2 + (\Delta_y - N v_{y0})^2 - N^2 (N+1)^2 \frac{A_{max}^2}{4} \le 0 \tag{5b}\label{NA5b}$$ $$(\Delta_x - N v_{x0})^2 + (\Delta_y - N v_{y0})^2 - N^4 \frac{A_{max}^2}{4} \le 0 \tag{5c}\label{NA5c}$$ depending on how the velocity is updated with respect to position updates, respectively. (This means each iteration when finding $N$ takes only something like 8 multiplications and 5 additions or subtractions. Since the number of iterations needed is roughly $2 \log_2 N$, this is actually pretty efficient. Even if we find $N \approx 10,000$, we only end up doing a couple of hundred multiplications and a hundred and thirty additions or subtractions. The algebraic solution has more terms than that!)
Furthermore, the optimized form is defined even for $N = 0$ and $N = 1$, so there is no risk of accidental divide-by-zero error. This makes the bisection search easier to write, too.
If we use
EQ5(N)for $\eqref{NA5a}$/$\eqref{NA5b}$/$\eqref{NA5c}$, the pseudocode to find $N$ is something likeWhen we do have $N$, it is just a matter of plugging it back to $\eqref{NA3a}$/$\eqref{NA3b}$/$\eqref{NA3c}$ to obtain the acceleration vector $(a_x , a_y)$.
The first While loop finds the smallest range containing a root, so that
Nmaxis a sufficient number of steps (yields an acceleration that does not exceed the given limit), andNminis an insufficient number of steps (yields an acceleration that exceeds the given limit). For each iteration $1 \le i \in \mathbb{N}$,Nmin= $2^{i-1}$ (exceptNmin=0for $i = 1$, andNmin=2for $i = 2$) andNmax=$2^i$, so only a few iterations are ever needed.The second While loop uses the bisection method (or binary search) to find the root within that range. As it halves the range on each iteration, it does at most as many iterations as the first loop did.
Because we want an $N$ that yields $a_x^2 + a_y^2 \le A_{max}^2$, we exclude the minimum (
Nmin) and include the maximum (Nmax) in the search. This is why we useN = Nmax - floor((Nmax - Nmin)/2), too: it never yieldsNmin. This is also why we start withNmin = 0: the smallestNreturned is 1.Here is an example Python (works in both Python 2 and Python 3) program, that solves the number of steps and the acceleration, when given the starting x and y coordinates, target x and y coordinates, initial velocity x and y components, and the maximum acceleration allowed:
The
solveacceleration()function takes the starting point, target or finish point, and the velocity as two-component tuples. If you include the commented outstderr.write()lines, the function will output a line per iteration when looking for the answer; you'll see that it does very few iterations even for the most complex cases. It returns the number of steps, plus the acceleration as a two-component tuple.The program itself outputs the parameters and the solution, including the location and velocity at each step, in a format suitable for e.g. gnuplot. If you save the output as e.g
out.txt, you can usein gnuplot to draw the trajectory, with small circles around each time step; hovering over each circle shows you the time step as a tooltip. To include the velocity in the tooltip, use
As it is written above, the example uses the $\eqref{NA2a}$ ("A") logic, position updated before velocity. I've commented the sections where you need to modify the code for the other two ("B" or "C").