Accelerating a point in two dimensions

Question

I have the following Python code:

import math

x = 0
y = 0

acceleration = 10
angle = 0

vx = 0
vy = 0

time = 10
for _ in range(time):
    vx += acceleration * math.cos(math.radians(angle))
    vy += -acceleration * math.sin(math.radians(angle))

    x += vx
    y += vy

print(x, y)

Which outputs:

550.0 0.0

This is not what the equation for displacement yields.

(acceleration * time**2) / 2 = 500

What am I doing wrong? I would like to solve the problem without using time; pretend it doesn't exist.

@Shubham Shaswat: Perhaps I should clarify: this is merely illustrative of the problem. In reality, I have an object which is updated every time step. — BathNinja, Feb 09 '20 at 14:59
The formula is exact, while you make a linear approximation. (Serious) rounding errors are inevitable. — dedObed, Feb 10 '20 at 14:33
`angle` never changes; you can precompute the sine and cosine outside of the loop. — chepner, Feb 10 '20 at 14:35
@dedObed: This is not a rounding error. This is a logical error. — BathNinja, Feb 11 '20 at 21:57
@chepner: A keen observation. Is it, however, useful to optimise useless code? — BathNinja, Feb 11 '20 at 21:59

dedObed · Accepted Answer · 2020-02-12T08:24:02.487

What you are trying to achieve is to find the exact integral of velocity over time, where velocity itself is given implicitly as the integral of acceleration. And you try do so by the simplest method available: the Euler method. Accumulating inaccuracies are inevitable.

On top of errors (imprecision) inherent to Euler method, your implementation has the error of updating variables in a sequential manner. i.e.: you combine past displacement with current velocity -- instead of with the corresponding past velocity. You should compute new values of each variable and update them at the same time. For example like this (omitting constants from your code):

import math                                                                                                                                                                                                        

acceleration = 10                                                                                                                                                                                                  
vx = 0                                                                                                                                                                                                             
x = 0                                                                                                                                                                                                              

for _ in range(10):                                                                                                                                                                                                
    new_x = x + vx                                                                                                                                                                                                 
    new_vx = vx + acceleration                                                                                                                                                                                     

    x = new_x                                                                                                                                                                                                      
    vx = new_vx                                                                                                                                                                                                    

print(x) # 450

In your current setup (with the fix), the simulation runs like this:

You can get better result by increasing the time resolution, e.g. by making steps of 0.1 instead of 1, you get:

If you're interested in better numerical integration methods, follow wikipedia to Runge-Kutta or Adams-Bashfort.

Here is the code to reproduce the plots:

import numpy as np                                                                                                                                                                                                 
import matplotlib.pyplot as plt                                                                                                                                                                                    

acceleration = 10                                                                                                                                                                                                  

t0 = 0                                                                                                                                                                                                             
t1 = 10                                                                                                                                                                                                            
nb_steps = 11                                                                                                                                                                                                      

ts = np.linspace(t0, t1, num=nb_steps)                                                                                                                                                                             
vs = np.zeros_like(ts)                                                                                                                                                                                             
xs = np.zeros_like(ts)                                                                                                                                                                                             

vs[0] = 0                                                                                                                                                                                                          
xs[0] = 0                                                                                                                                                                                                          

true_xs = acceleration * ts ** 2 / 2                                                                                                                                                                               


for i, t in enumerate(ts):                                                                                                                                                                                         
    if i == 0:                                                                                                                                                                                                     
        continue # initial conditions are preset                                                                                                                                                                   

    delta_t = t - ts[i-1]                                                                                                                                                                                          
    vs[i] = vs[i-1] + acceleration * delta_t                                                                                                                                                                       
    xs[i] = xs[i-1] + vs[i-1] * delta_t                                                                                                                                                                            

plt.figure()                                                                                                                                                                                                       
plt.plot(ts, vs, label='velocity')                                                                                                                                                                                 
plt.plot(ts, xs, label='displacement-sim')                                                                                                                                                                         
plt.plot(ts, true_xs, label='displacement-true')                                                                                                                                                                   
plt.legend()                                                                                                                                                                                                       
plt.show()

The semi-implicit or symplectic Euler method is also an order 1 method and in some situations like orbital mechanics qualitatively more correct than the explicit or implicit Euler methods. It is only a half of a half-step away from the Verlet method, which is even better for mechanical systems. So it is right to say that the question code failed to implement the explicit Euler method, but wrong to say that it is categorically wrong. — Lutz Lehmann, Feb 13 '20 at 00:00
Thank you. This is an excellent answer. I had no idea there were methods to these things; I assumed they were precise. I have only one further thing to ask: In my code, when I set starting velocity to `vx = -acceleration * math.cos(math.radians(angle)) / 2` and `vy = acceleration * math.sin(math.radians(angle)) / 2`, I seem to get correct and precise (x, y) results (at the cost of incorrect velocity). Why is that? — BathNinja, Feb 13 '20 at 15:52

Shubham Shaswat · Answer 2 · 2020-02-10T14:30:39.200

0

In your case the x and y should be updated with average of vx initial and final velocity and vy respectively.

As the

If you want x and y for every time step you should do like this,

import math

x = 0
y = 0

acceleration = 10
angle = 0

vx = 0
vy = 0

time = 10

vx = (vx + acceleration * math.cos(math.radians(angle))*time)/2 #average of velocity
vy = (vy  -acceleration * math.sin(math.radians(angle))*time)/2 #average of velocity

for _ in range(time):

  x += vx
  y += vy

  print(x, y)

edited Feb 10 '20 at 14:30

answered Feb 09 '20 at 15:04

Shubham Shaswat

1,250
9
14

I have updated my request, in that I do not want to use 'time' to solve the problem. You did not answer the question "What am I doing wrong?" You just wrapped the equation for displacement in a loop. – BathNinja Feb 09 '20 at 17:46

Accelerating a point in two dimensions

2 Answers2