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.

Answer 1

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()