# Calculate the constant acceleration needed for the discrete time trajectory to intersect a given target point

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?

• It looks like a variational problem. Any conditions over the path to the goal? Minimum distance? – Cesareo Aug 11 '18 at 13:35
• When you write “constant acceleration,” do you also mean for the direction of the acceleration to be constant? Does the time at which the object reaches its target location matter? – amd Aug 11 '18 at 17:27
• Also, your formula for the total distance traveled is really only valid for straight-line motion. If the object is changing direction, then some part of that total acceleration is applied in a direction perpendicular to the path, which reduces the amount of acceleration that can change the object’s speed, and so also the distance it travels. For an extreme example, take uniform circular motion, in there’s a constant amount of acceleration that’s always perpendicular to the (instantaneous) direction of travel: this acceleration has no effect whatsoever on the object’s speed. – amd Aug 11 '18 at 17:29
• @Cesareo No, only that the amount added to the velocity in each step is the same and does not exceed the maximum acceleration. – Misys Aug 11 '18 at 21:01
• If the magnitude and direction of the acceleration is fixed, we have a case of projectile motion. The "ground" would be a line passing through the target point. The direction of this line would depend on the direction of the initial velocity. The "horizontal" component of this velocity would need to be toward the target point. A useful formula combining angle, acceleration and target position is found here. – Jens Aug 12 '18 at 2:02

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 like

Function StepsNeeded(x0,y0, xN,yN, vx,vy, Amax):
Let  dx = xN - x0
Let  dy = yN - y0
Let  Nmin = 0
Let  Nmax = 2

# Find the Nmin,Nmax range spanning the solution
While EQ5(Nmax) > 0:
Let  Nmin = Nmax
Let  Nmax = Nmax * 2
End While

Let  N = Nmax
While Nmax > Nmin:
Let  N = Nmax - floor((Nmax - Nmin)/2)
Let  C = EQ5(N)
If C > 0:
Let  Nmin = N
Else If C < 0:
If Nmax > N:
Let  Nmax = N
Else:
Break While
End If
Else:
Break loop
End If
End While

Return N
End Function


When 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 Nmax is a sufficient number of steps (yields an acceleration that does not exceed the given limit), and Nmin is 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}$ (except Nmin=0 for $i = 1$, and Nmin=2 for $i = 2$) and Nmax=$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 use N = Nmax - floor((Nmax - Nmin)/2), too: it never yields Nmin. This is also why we start with Nmin = 0: the smallest N returned 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:

from math import sqrt
from sys import stderr, argv, exit

# This file is in public domain. No guarantees, no warranties.
# Written by Nominal Animal <question@nominal-animal.net>.

def solveacceleration(start, finish, velocity, maxaccel):
dx = finish[0] - start[0]
dy = finish[1] - start[1]
vx = velocity[0]
vy = velocity[1]
ac = 0.25 * maxaccel * maxaccel

if maxaccel <= 0.0:
return 0, (0.0, 0.0)

nmin = 0
nmax = 2
while True:
cx = dx - nmax*vx
cy = dy - nmax*vy

# A) cc = cx**2 + cy**2 - ac * nmax**2 * (nmax-1)**2
# B) cc = cx**2 + cy**2 - ac * nmax**2 * (nmax+1)**2
# C) cc = cx**2 + cy**2 - ac * nmax**4

cc = cx**2 + cy**2 - ac * nmax**2 * (nmax-1)**2
# stderr.write("nmin=%u, nmax=%u, cc=%.6f\n" % (nmin, nmax, cc))
if cc > 0:
nmin = nmax
nmax = nmax * 2
else:
break

n = nmax
while nmax > nmin:
n = nmax - int((nmax - nmin) / 2)
cx = dx - n*vx
cy = dy - n*vy

# A) cc = cx**2 + cy**2 - ac * n**2 * (n-1)**2
# B) cc = cx**2 + cy**2 - ac * n**2 * (n+1)**2
# C) cc = cx**2 + cy**2 - ac * n**4
cc = cx**2 + cy**2 - ac * n**2 * (n-1)**2

# stderr.write("nmin=%u, nmax=%u, n=%u, cc=%.6f\n" % (nmin, nmax, n, cc))
if cc > 0:
nmin = n
elif cc < 0:
if nmax == n:
break
else:
nmax = n
else:
break

# A) a_ = 2*(d_ - n*v_) / (n * (n - 1))
# B) a_ = 2*(d_ - n*v_) / (n * (n + 1))
# C) a_ = 2*(d_ - n*v_) / (n**2)
ax = 2*(dx - n*vx) / (n * (n - 1))
ay = 2*(dy - n*vy) / (n * (n - 1))
return n, (ax, ay)

if __name__ == '__main__':

if len(argv) != 8:
stderr.write("\n")
stderr.write("Usage: python %s x0 y0 x1 y1 xv yv amax\n" % argv[0])
stderr.write("\n")
exit(0)

p0 = ( float(argv[1]), float(argv[2]) )
p1 = ( float(argv[3]), float(argv[4]) )
v0 = ( float(argv[5]), float(argv[6]) )
amax = float(argv[7])

n, a = solveacceleration(p0, p1, v0, amax)

print("#  From: %9.3f %9.3f" % p0)
print("#    To: %9.3f %9.3f" % p1)
print("#    v0: %+9.3f %+9.3f" % v0)
print("# Steps: %9d" % n)
print("#     a: %+9.3f %+9.3f  (%.3f of maximum %.3f)" %
(a[0], a[1], sqrt(a[0]*a[0]+a[1]*a[1]), amax))

p = p0
v = v0
print("# step    x         y          vx        vy")
print("")
for i in range(0, n+1):
# A) print("%-5d %9.3f %9.3f   %+9.3f %+9.3f" % (i, p[0], p[1], v[0], v[1]))
#    p = (p[0]+v[0], p[1]+v[1])
#    v = (v[0]+a[0], v[1]+a[1])

# B) v = (v[0]+a[0], v[1]+a[1])
#    print("%-5d %9.3f %9.3f   %+9.3f %+9.3f" % (i, p[0], p[1], v[0], v[1]))
#    p = (p[0]+v[0], p[1]+v[1])

# C) v = (v[0]+0.5*a[0], v[1]+0.5*a[1])
#    print("%-5d %9.3f %9.3f   %+9.3f %+9.3f" % (i, p[0], p[1], v[0], v[1]))
#    p = (p[0]+v[0], p[1]+v[1])
#    v = (v[0]+0.5*a[0], v[1]+0.5*a[1])

print("%-5d %9.3f %9.3f   %+9.3f %+9.3f" % (i, p[0], p[1], v[0], v[1]))
p = (p[0]+v[0], p[1]+v[1])
v = (v[0]+a[0], v[1]+a[1])


The solveacceleration() function takes the starting point, target or finish point, and the velocity as two-component tuples. If you include the commented out stderr.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 use

plot "out.txt" u 2:3 notitle w lines lc -1, \
"out.txt" u 2:3:1 notitle w labels \
hypertext point pt 6


in 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

plot "out.txt" u 2:3 notitle w lines lc -1, \
"out.txt" u 2:3:(sprintf("%s: (%s,%s)", stringcolumn(1), stringcolumn(4), stringcolumn(5))) \
notitle w labels hypertext point pt 6


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").

• Thanks! I was trying to avoid having to search for the solution, but I didn't realize how complicated the algebraic solutions would be. – Misys Aug 14 '18 at 6:13
• @Misys: I do have code implementing the discrete time solution in C. Let me know if you'd like me to post a C or Python example that finds $N$ (and therefore also $a_x$ and $a_y$), given starting point $(x_0, y_0)$, target $(x_N , y_N)$, initial velocity $(v_x , v_y)$, and maximum allowed acceleration $A$ (or squared, $A^2$). – Nominal Animal Aug 14 '18 at 19:04
• Sure! No need for me to reinvent the wheel if you already have one. I'm more familiar with python, but if you already have it written in C that works too. – Misys Aug 15 '18 at 6:56
• @Misys: Added; the code is Public Domain, so use as you wish, but don't blame me if it does not work. (I did test it, lightly.) If you care to modify your question slightly, to mention discrete time (say, "Calculate the constant acceleration needed for the discrete time trajectory to intersect a given target point"), I could simplify my answer quite a bit, and omit the continuous-time initial part completely. (Where continuous time trajectories use differential equations and derivatives, discrete time uses difference equations.) – Nominal Animal Aug 16 '18 at 0:49
• The reason for that suggestion is, that I believe that way others working on a similar problem would be more likely to encounter your question, and we could help those future people too, by presenting the matter with a slightly better focus. I think. ;) – Nominal Animal Aug 16 '18 at 0:53

Your formula for distance traveled under acceleration works if the acceleration is added to the velocity at the beginning of each time step. If the acceleration is added at the end of each step, you should change $n(n+1)$ to $n(n-1).$ But if your physics engine is cleverly implemented, it might “split the difference” in effect, in which case the correct formula would have $n^2$ instead of $n(n+1).$

Any of these possible implementations is roughly the same difficulty to solve. Consider two vectors, one representing the object’s total travel after $n$ steps assuming no acceleration, and the other representing the total travel under your chosen acceleration assuming no initial velocity. The object’s actual travel is the sum of these vectors. You want that sum to equal the vector from the object’s initial position to the target.

Let $\theta$ be the angle between the initial velocity and the initial direction to the target. The vector sum forms a triangle with sides $d$ (the initial distance to the target), $vn$ (in the initial direction of travel), and $\frac a2 n(n+1)$ opposite the angle $\theta.$ Use the Law of Cosines: $$\left(\frac a2 n(n+1)\right)^2 = v^2 n^2 + d^2 + 2dvn \cos\theta.$$ This comes out to a fourth-degree polynomial in $n,$ which has an “exact” solution in principle but is usually solved by numeric methods (basically various more or less sophisticated forms of “guess and check”).

But if you will settle for a solution that is not necessarily optimal, you can solve the following equation instead: \begin{align} \left(\frac a2 n(n+1)\right)^2 &= v^2 n^2 + d^2 + 2dvn \\ &= (vn + d)^2. \end{align} This is easily solved by taking the square root of both sides (knowing they both must be positive) and then solving a quadratic equation.

Since $2dvn \geq 2dvn \cos \theta,$ the simplified equation will give you a number of timesteps $n$ that is certainly sufficient, and possibly more than you want at maximum acceleration. In the latter case (that is, whenever $\theta\neq0$), after solving the simplified equation for $n,$ plug that value of $n$ into the original equation and replace $a$ by a factor that makes both sides equal. Finally, use that factor as your acceleration. You can find the correct direction of acceleration by recomputing the vector sum and finding the direction of the vector representing acceleration.