Construct an Euler’s method code for simple one-dimensional motion.
Background Information
Review of Kinematics
Most will be familiar with the kinematic equations of motion, which can be used in situations involving constant acceleration. The kinematic equations are. \[v_f = v_i + a \Delta t\]\[x_f = x_i + v_i \Delta t + {1\over 2} a \Delta t^2\]
As an example, consider a car that is initially at rest, waiting for the stoplight to turn green. Once the light turns green, the driver presses on the accelerator and the car’s speed begins to increase. After accelerating at a rate of \(a = 5\) m/s\(^2\) for \(t = 3\) s, the driver eases up and reduces the acceleration to a more modest \(a = 2\) m/s\(^2\) and maintains this acceleration for another \(\Delta t = 8\) s. How would you determine the position and speed of the car after the full \(11\) seconds had elapsed? The kinematic equations above could be applied to the first 3 seconds of the trip or to the last 8 seconds of the trip, but couldn’t be applied from the start to the end. Let’s use the kinematic equations to determine the speed and position after the first 3 seconds. \[ v_1 = v_0 + a_0 \Delta t\]\[ = (0 ~\text{m/s}) + (5 ~\text{m/s}^2) (3 ~\text{s})\]\[ = 15 ~\text{m/s}\]\[x_1 = x_0 + v_0 \Delta t + {1\over 2} a_0 \Delta t^2\]\[ = (0 ~\text{m}) + (0 ~\text{m/s}) (3 ~\text{s}) + {1\over 2}(5 ~\text{m/s}^2) (3 ~\text{s})^2\]\[ = 22.5 ~\text{m}\] Now that we have the speed and velocity at \(t = 3\) we can proceed to use the same kinematic equations for the second part of the motion. \[ v_2 = v_1 + a_1 \Delta t\]\[ = (15 ~\text{m/s}) + (2 ~\text{m/s}^2) (8 ~\text{s})\]\[ = 31 ~\text{m/s}\]\[x_2 = x_1 + v_1 \Delta t + {1\over 2} a_1 \Delta t^2\]\[ = (22.5 ~\text{m}) + (15 ~\text{m/s}) (8 ~\text{s}) + {1\over 2}(2 ~\text{m/s}^2) (8 ~\text{s})^2\]\[ = 206.5 ~\text{m}\] Since the kinematic equations only apply to constant-acceleration motion, the problem had to be divided into two parts: one for each constant-acceleration portion of the motion.
Now, what if the motion involved 5 different segments, each with a different acceleration, with the acceleration starting at 4 m/s\(^2\) and decreasing by \(1\) m/s\(^2\) every second.
time(s)
acceleration (m/s\(^2\))
0
4
1
3
2
2
3
1
We could repeat the process from above to generate the positions and velocities at these four times.
The \(x_{n+1}\) and \(x_n\) are standard notation to mean the next x (\(x_{n+1}\)) and the current x (\(x_n\))
This process could be extended for as long as you wanted but I hope you can see that you don’t want to do this by hand for very long.
Euler’s Method
Most realistic scenarios involve accelerations that change continuously, not discretely like the examples above. However, we can still employ the method used above if we use a time step (\(\Delta t\)) that is very small. Even if the acceleration varies substantially over long time intervals, it will always be nearly constant over very small time steps. We call this Euler’s method and the steps to implement Euler’s method are as follows:
Define the initial conditions. For particle motion this involves defining the initial position and velocity of the object as well as any other variable that do not change over time.
Repeatedly calculate the following quantities.
The current acceleration. The acceleration function is typically given to you or can be derived from Newton’s law. This expression will always be a function of time,position, and/or velocity and must be re-evaluated every time you advance by a time step.
Use a kinematic equation to calculate the next position using the previous position and velocity.
Use a kinematic equation to calculate the next velocity using the previous position and velocity.
If you want to store all of your positions and velocities for later use, append them to a list.
An example implementation of Euler’s method for a particle (experiencing no air resistance) thrown upward (with \(v = 30\) m/s) from an initial height of \(y_i = 70\) m, is shown below.
# Import statementsimport numpy as np import matplotlib.pyplot as pltv0 =30.0# Initial velocity in m/sy0 =70# Initial height in mdelta_t =0.01# Time step in seconds (should be small)t0 =0# Start time in secondst_max =7.0#Final TimeN = (t_max - t0)/delta_t # Calculate the number of time stepsN =int(N) # Make sure N is an integer# Make lists to store your positions, velocities and times.y = [y0]v = [v0]t = [t0]for i inrange(N): a =-9.8# Calculate the current acceleration# Store the new position, velocity, and time in the lists v.append(v[i] + a * delta_t) y.append(y[i] + v[i] * delta_t) t.append(t[i] + delta_t)plt.plot(t,y,linewidth =2, color ='red')plt.xlabel("time (s)")plt.ylabel("Position (m)")plt.title("Position vs. time for a ball in free fall")
Text(0.5, 1.0, 'Position vs. time for a ball in free fall')
The code - Broken down
This section will break the program down piece by piece. Here’s the first section:
"""One Dimensional free-fall Euler CodePH150Brother NelsonThis code will calculate the exact solution for a ball in free fallhaving been shot straight up using Euler's method."""#Import numerical and plotting packagesimport numpy as npimport matplotlib.pyplot as plt
Everything inside the triple quotes (“““) counts as one long comment. If you ever need to add comments that will take more than one line, you can start and end it with """. It’s a great way to leave yourself notes about what each program does. Below the introductory comment, numpy and matplotlib are imported to do the number crunching and plotting later on.
This next part sets up all of the constants and inputs. It’s a good idea to put it at the beginning of the program so that they are easy to change.
#Initial conditions and physical setup constantsv0=30.0#Initial velocity in m/sx0=70#Initial Postion in m#Set up the time steps and number of calculationsdelta_t =0.01#Time step in secondst0=0#Start time in secondst_max =7.0#Final time
The part with x0 and v0 are the starting position and velocity. delta_t is time step or amount of time that will advance for each calculation. Remember, Euler’s method works by assuming that acceleration is roughly constant over short time intervals. delta_t sets that time interval.
With initial values defined, the number of calculations that will be needed can be determined. Here’s the piece of code that does that:
#Calculate the number of timesteps we need to makeN=(t_max-t0)/delta_t#Make sure N is an integerN=int(N)
Python keeps two types of numbers: Floats - numbers with decimal points and Integers - whole numbers. For example: \(2.0\) is a float, \(2\) is an integer. Even though they have the same value, Python sees them as different things. You can only get an item from a list using an integer. (myList[2] works, but myList[2.0] won’t). The command N=int(N) takes whatever number N is and converts it to an integer. The int command truncates instead of rounding, so int(2.9) becomes 2, not 3.
Next, lists that will eventually store all of the positions, velocities, and times are initialized. Putting square brackets around something tells Python that it is a list. Since more x,v, and t values are going to be added later, we need to warn Python that these are going to be lists. Initially, these lists only contain the initial values.
#Make lists to hold our positions, velocities, and timesx=[x0]v=[v0]t=[t0]
This next section of code is what is called a loop:
#Now perform an Euler's method calculation.for i inrange(N):#Find the current acceleration a=-9.8# Store the new position, velocity, and time in the lists v.append(v[i] + a * delta_t) y.append(y[i] + v[i] * delta_t) t.append(t[i] + delta_t)#End of Euler loop =========================================
Breaking down the for loop
For most of you, a loop will be a new idea. Loops are very helpful when you need to tell a computer to do something over and over again. Here’s the general structure:
for<Thing that Changes>: Stuff the computer does to/with the thing that changesNot part of the loop
The code that gets inserted at: <Thing that Changes> should be a list of items that are to be iterated over. As an example, consider a famous mathematical sequence of numbers called the Fibonacci sequence. The first few numbers in the Fibonacci sequence are: [1,1,2,3,5,8,13] (Can you see the pattern). A loop to iterate over this list of numbers is shown below:
fibonacci = [1,1,2,3,5,8,13]for num in fibonacci:print(num)
1
1
2
3
5
8
13
When this code is executed, the sequence of numbers first gets stored in the variable fibonacci. The for num in fibonacci tells Python to look at each thing in the list one by one, and call it num, each in turn. Therefore, this program would first print 1, then 1 (again), then 2, etc. Once all of the numbers have been printed one at a time, the final command will execute and print the entire list. This last statement won’t execute until the loop has successfully iterated over everything in the list named fibonacci.
Now let’s look at the for statement in the Euler’s code: i in range(N). The range(N) command makes a list of integers from zero to N-1. If N was five, range(N) would be [0,1,2,3,4]. Therefore, the first time through the for loop i=0, the second time i=1, and so on until we reach N-1.
Now let’s look at the inside of the loop:
#Find the current acceleration a=-9.8# Store the new position, velocity, and time in the lists v.append(v[i] + a * delta_t) y.append(y[i] + v[i] * delta_t) t.append(t[i] + delta_t)
This is where Euler’s method is found. First the acceleration is calculated, which is a constant -9.8 in this case but usually will be a function of time, position, and/or velocity. The next statement involves an append statement:
v.append(v[i] + a * delta_t)
The append statement is used to add something to the end of a list. In this case, the statement v.append(v[i]+a*delta_t) adds v[i]+a*delta_t to the end of the list named v. This is nothing more than the kinematic equation \(v_f=v_i+a\Delta t\), using v[i] as \(v_i\) to calculate the “final” velocity. Each time the computer goes through the loop, a different value for i is used, and therefore a different value for v[i]. If the list starts with just one value in it (v=[0]), then by the end of the first iteration there are two things in the list. (The same can be said for the position and time lists)
The first few iterations of this loop are done by hand and placed in the table below to help you see exactly what the loop is doing.
time
$$y_f = y_i + v_i \Delta t$$
$$v_f = v_i + a \Delta t$$
0
0.0
70.000
30.00
1
0.1
73.000
29.02
2
0.2
75.902
28.04
3
0.3
78.706
27.06
Velocity can be used to calculate the position using the kinematic equation \(x_f=x_i+v_i\Delta t+\frac{1}{2}a\Delta t^2\). The term \(\frac{1}{2}a\Delta t^2\) is neglected because if \(\Delta t\) is small then \(\Delta t^2\) will be very small, making this entire term inconsequential. Once the loop is finished, the program just builds a plot of the data:
#Plot the resultsplt.plot(t,x,linewidth=2,color='red')plt.xlabel('time (s)')plt.ylabel('Position (m)')plt.title('Position vs. time for a ball in free fall')
Air Drag
This week we will be using Euler’s method to predict the fall time for an object experiencing a non-negligible amount of drag. The drag force acting on a falling object is given by:
\[F_D = {1\over 2} \; C \rho A v^2 \;\;\;\;\;(1) \]
where \(C\) is the drag coefficeint and has the following values: 0.5 for a sphere and 0.3 for a cone.
\(\rho\) is the air density and has a value of 1.02 kg/m\(^3\) in Rexburg
A is the cross-sectional area
v is the velocity of the object
Applying Newton’s second law to the free-body diagram shown above gives:
Solving for the acceleration in the y-direction gives:
\[ a_y = {d v_y \over dt} = \frac{ C \rho A v^{2}}{2m} - g \;\;\;\;\;(3)\]
Notice that this expression for the acceleration is not constant; it changes as \(v\) changes. The kinematic equations can be used to update the position and velocity \[ y_{n+1} = y_{n} + v_{n} \cdot \Delta t \;\;\;\;\;(5)\]
\[ v_{n+1} = v_{n} + a_{n} \cdot \Delta t \;\;\;\;\;(6)\]
\[ = v_{n} + ({C \rho A v^2 \over 2 m} - g) \cdot \Delta t \;\;\;\;\;(6)\]
Python skills
In this lab you will need to use programming loops to use Euler’s method. Loops allow programs to rerun the same block of code multiple times. This might seem like a funny thing to want to do but it turns out that there are many important tasks that are repetitive in nature (perhaps with small changes for each successive repetition). A loop provides a succinct and efficient way to perform tasks of this nature.
For loop
The for loop is probably the most common loop you will encounter and is a good choice when you know beforehand exactly what things you want to loop over. Here is an example of a for loop that is used to add up the elements of a list.
which isn’t that much longer than using a loop. However, as the list gets longer and/or the mathematical operations being performed get more complex the second method would get unreasonably long.
The correct language is to say that we are iterating over the list [3,2,1,9.9]. This means that the loop variable (i in this case but you can choose it to be whatever you want) gets assigned the values of the list elements, one by one, until it reaches the end of the list. You can use for loops to iterate over any multi-element object like lists or tuples. Python uses indentation to indicate where the loop ends. In this case there was only one statement inside to loop, but if you wanted more than one each line should be indented.
Often you will want to iterate over a list of integers. The range function is a good choice for this. The range function takes up to 3 arguments. The first is the starting number, the second is the ending number, and the third is the step size.
range(start,stop,stepsize)
for i inrange(5,50,3): #Generates a list from 5 -> 50 with a step size of 3print(i)
5
8
11
14
17
20
23
26
29
32
35
38
41
44
47
These examples are so simple that you might wonder when a loop might actually be useful to you. Let’s see if we can build a loop to calculate the following sum
\[ \sum_{n=1}^{1000} {1\over n^2} \tag{1}\]
theSum =0for n inrange(1,1000): theSum = theSum +1/n**2print(theSum)
1.6439335666815615
Here, n is being assigned the values 1,2,3,4....1000, one by one, until it gets all the way to 1000. Each time through the loop, n is different and the expression 1/n**2 evaluates to a new value. The variable theSum is updated each time through to be the running total of all calculations performed thus far. Here’s another example of a loop used to calculate the value of \(20!\):
theProduct =1for n inrange(1,21): theProduct = theProduct * n #Multiply theProduct by nprint(theProduct)
2432902008176640000
Remember that the range function creates a list starting at \(1\), going up to \(21\) but not including it. The math library has a function called factorial that does the same thing. Let’s use it to check our answer:
from math import factorialfactorial(20)
2432902008176640000
Boolean logic inside loops
Often when using loops, we only want a block of code to execute when some condition is satisfied. We can use boolean logic inside of the loop to accomplish this. For example, let’s write a loop to compute the following sum:
\[ \sum_{{n\over 5} \in \text{ Int and } {n\over 3} \in \text{ Int}} {1\over n^2} \]
which is similar to the one we did above, but this time we only want to include terms where \(n\) is a perfect multiple of both 5 and 3. To check to see if n is a perfect multiple of a number we can calculate the modulo (remainder after division) using the % operator and check that it is equal to zero.
theSum =0for n inrange(1,1000):if n %5==0and n %3==0: theSum = theSum +1/n**2print(theSum)
0.007243985583159138
Activity I: Minimum Fall Time
Goal (Overview)
Using kinematics, calculate the time it takes for a particle to fall, starting from rest, through a distance of \(6\) meters, neglecting air drag. This number will serve as a point of comparison when we include air drag in our analysis.
Procedure
The kinematic equations are given below. \[y_f = y_i + v_i \Delta t + {1\over 2}a \Delta t^2\]\[v_f = v_i + a \Delta t\]\[v_f^2 = v_i^2 + 2 a \Delta x\] Use these equations to calculate the time it takes for a particle to fall through a distance of \(6.0\) meters starting from rest. Neglect air resistance. Show your math and calculations below.
Include kinematic details here
# Code to perform calculations here
Activity II: Modeling the motion using Euler’s method
Goal (Overview)
Write some python code that will accurately predict the fall time for a projectile that is experiencing air resistance.
Procedure
In the example discussed above, a model for a particle that is thrown upward and experiences no air resistance was constructed. Use equations 3-6 from the “Air Drag” section to modify that code to include air resistance. Set \(\rho = 1.02\) kg/m\(^3\), \(C = 0.5\) and make reasonable guesses for the values of \(A\) and \(m\). We’ll fine tune these values later.
Modify the code to model a particle that is dropped from rest rather than thrown upward initially.
Your value for \(\Delta t\) will need to be reasonably small to get accurate results. Decrease your value of \(\Delta t\) until the results stop changing.
To find the flight time, you should use trial and error on t_max until the particle hits the ground.
As a test of the correctness of your code, increase the initial height of the projectile and verify that the drop time changes as you would expect.
As a test of the correctness of your code, increase the mass of the projectile and verify that the drop time changes as you would expect.
As a test of the correctness of your code, decrease the drag constant of the projectile and verify that the drop time changes as you would expect.
# Import from librariesfrom numpy import pi# Define constants# Set initial conditions# Loop to calculate flight time using Euler's methodwhile# Put your stopping condition here:plt.ylim(0,y0)
Activity III: Predicting times
Procedure
Make mass and diameter measurements on each of the objects given to you.
Choose a location to perform the drop and measure the initial height. You should find a location that provides a large drop distance. (The track above the I-center or the foyer in the Manwaring center are good choices. We’ve dropped them from the top of the stadium before but a tall fence was recently installed which makes it hard to drop from there. )
Then use the code you wrote above to predict the drop times for each of the six objects.
Fill in the table below with your predictions.
Object
Mass (kg)
Diameter (m)
Precited fall time (s)
Ping Pong Ball
Styrofoam cone
Styrofoam sphere (small)
Styrofoam sphere (medium)
Beach ball
Activity IV: Comparison of Predicted and Measured Times
Measure the flight time for each object \(10\) times.
