Guide to Euler-Cromer Method¶
Those of you who had Matter and Interactions for your Intro Physics book will be mostly familiar with this method. Nothing really changes here from how you did things previously, I'll just show a little more nuanced ways to doing the calculation. Euler-Cromer can be particularly useful in situations where you might have different forces switching on or off throughout the motion that would be difficult to capture using Scipy's Odeint.
The Math¶
The idea behind the Euler method is to use small time steps to approximate the infinitesimally small time steps of a derivative. Then we slowly iterate through a solution, plugging the results from the previous computation into the next computation. Starting with Newton's 2nd law: \begin{align} a &= \frac{F}{m} \\ \dot v &= \frac{F}{m}\\ \end{align} We then approximate this as \begin{align} \frac{\Delta v}{\Delta t} &= \frac{F}{m} \\ \Delta v &= \frac{F}{m}\Delta t \\ v_f &= v_i + \frac{F}{m}\Delta t \tag{$\dagger$} \end{align} And thus we have a method to predict a future velocity after a small time increment $\Delta t$. We can repeat this to get back to the position: \begin{align} \dot r &= v \\ \frac{\Delta r}{\Delta t} &= v \\ \Delta r &= v\Delta t \\ r_f &= r_i + v\Delta t \tag{$\ddagger$} \end{align} The key assumption (and the Cromer part of Euler-Cromer) is that the instantaneous velocity in ($\ddagger$) is approximately equal to the final velocity in ($\dagger$). This is pretty reasonable as long as the time step is small enough that the velocity does not change that much over the time interval! So when in doubt, err on the side of making your time intervals too small when using this method!
Basic Usage¶
Regardless of how you do Euler-Cromer, you'll need to loop through the solution. Whether you want to do things with a while loop or a for loop is mostly up to you, but I'll show you both versions here. First, consider the problem:
Imagine the case of a projectile flying through the air, under the influence of a constant gravitational force. Say we launch the 10 kg projectile with an initial velocity of 50 meters per second at an angle of $60^\circ$ above the horizontal.
We'll start by declaring our constants and importing numpy. dt here is the duration of each time step, so I'm taking a pretty large 0.1 second per step in this case.
import numpy as np
m = 10 # kg
g = 9.8 # m/s^2
vi = 50 # m/s
a = 60*np.pi/180 # radians
t = 0 # second
dt = 0.1 # second
Using while loops:¶
Frequently we'll want to compute the solution up to a certain threshold. For example, for 10 seconds or until the object returns to the ground. This type of condition lends itself nicely to while loops. We could compute the desired trajectory of our projectile then using the below code:
# Initial conditions
vx = vi*np.cos(a)
vy = vi*np.sin(a)
x = 0
y = 0
while t < 0.5:
# Calculate the force
F = -m*g
# Calculate the acceleration
ax = 0
ay = F/m
# Update the velocity
vx = vx + ax*dt
vy = vy + ay*dt
# Update the position
x = x + vx*dt
y = y + vy*dt
# Update the time
t = t + dt
print(x,y)
And you see we get a printout of the x position and y position for each bit of time. Unfortunately, we are overriding these values every iteration. This is fine if you just want the end result, but if you want the results throughout the flight, for example for graphing purposes, we need a way to retain the data while still iterating forwards. The best way I know to do this is by appending each iterations results to a growing list. This is simple to implement:
m = 10 # kg
g = 9.8 # m/s^2
vi = 50 # m/s
a = 60*np.pi/180 # radians
t = 0 # second
dt = 0.01 # second
# Initial conditions
vx = vi*np.cos(a)
vy = vi*np.sin(a)
x = 0
y = 0
# Initial lists of empty x and y values
xs = []
ys = []
while t < 10:
# Calculate the force
F = -m*g
# Calculate the acceleration
ax = 0
ay = F/m
# Update the velocity
vx = vx + ax*dt
vy = vy + ay*dt
# Update the position
x = x + vx*dt
y = y + vy*dt
# Append these latest results to our lists
xs.append(x)
ys.append(y)
# Update the time
t = t + dt
Now, at the end, we have an actual list of all the $x$ and $y$ values that the object possessed throughout its flight, as can be seen by looking at the first 10 values of xs:
xs[:10]
Using for loops¶
If you want your solution over a particular chunk of time, in often times might make more sense to just prepopulate some of the lists early and then loop over them. This also cleans up some of the syntax (in my opinion) a little. To repeat the above, we could do something like:
# Constants
m = 10 # kg
g = 9.8 # m/s^2
vi = 50 # m/s
a = 60*np.pi/180 # radians
t = 0 # second
dt = 0.01 # second
# Generate all desired times
ts = np.arange(0,10,dt)
# For each time we'll need position and velocity in x and y directions
# We can just set these to 0 initially
xs2 = np.zeros(ts.shape)
ys2 = np.zeros(ts.shape)
vxs2 = np.zeros(ts.shape)
vys2 = np.zeros(ts.shape)
# Setting your initial conditions is the same as setting the first value of each array
xs2[0] = 0
ys2[0] = 0
vxs2[0] = vi*np.cos(a)
vys2[0] = vi*np.sin(a)
# And then we for loop our way through everything, implementing Euler-Cromer
for i in range(0, len(ts)-1):
# Calc Force
Fx = 0
Fy = -m*g
# Calc next velocity
vxs2[i+1] = vxs2[i] + Fx/m*dt
vys2[i+1] = vys2[i] + Fy/m*dt
# Calc next position
xs2[i+1] = xs2[i] + vxs2[i+1]*dt
ys2[i+1] = ys2[i] + vys2[i+1]*dt
Plot comparison¶
Just to show both solutions got us the same solution, we can check out a plot of the results:
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-darkgrid')
fig = plt.figure(figsize=(10,7))
plt.plot(xs,ys, lw=3, label='While loop solution')
plt.plot(xs2, ys2, '--', lw=3, label='For loop solution')
plt.xlabel('Horizontal distance (m)')
plt.ylabel('Vertical height (m)')
plt.legend()
plt.show()
