Before we start, here is a very simple tutorial on finite differences.
Heun's method is not explained in the textbook, and is defined on Wikipedia as the 'improved or modified Euler's method, or a similar two-stage Runge-Kutta method. [..] it is a numerical procedure for solving ODEs with a given initial value. Both variants can be seen as extensions of the Euler method into two-stage second-order Runge-Kutta methods.'
$h$ is step size. We're going to expand $f$ as a function of $h$ on the righthand side:
$$y(x + h) = y(x) + \frac{h}{2} f (x, y(x)) + \frac{h}{2} f (x + h, y(x) + hf (x, y(x)))$$Next, expanding the $h$ on the righthand side:
$$y(x + h) = y(x) + \frac{h}{2} \left( f (x, y(x)) + h \left( \frac{\partial f}{\partial x} + f (x, y(x)) \frac{\partial f}{\partial y} ) \right) + O(h^2) \right) $$This is an explanation on what the big O is. It is part of the Bachmann-Landau notations. This is a mathematical notation that describes the limiting behavior of a function when the argument tends towards a particular value or infinity.
\mathcal{O} gives $\mathcal{O}$
This corresponds with the given Taylor expansion in the textbook. It is a third-order approximation, in accordance with the midpoint method. A video explainer on Runge-Kutta can be found here.
In the book, it was given that: $$\frac{dy}{dx}=f(x,y)$$ $$\frac{dy}{dt}=-y$$
The Euler method was given:
$$y(x+2h)=y(x+h)+hf(x+h,y(x+h)$$Complementary explanation on Euler's method on Freecodecamp.
Evaluate how the average local error error in the final value and slope depend on step size.
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = 12, 6
def euler(y0, dy0, step, max_step):
y = y0
dy = dy0
result = [y]
t = step
while t <= max_step:
tmp = y
y = y + step * dy
dy = dy - step * tmp
result.append(y)
t += step
return np.array(result)
y = euler(1, 0, 0.001, 100 * np.pi)
plt.plot(y, color="pink")
plt.show()
Evaluate how the average local error error in the final value and slope depend on step size
The point of using an adaptive stepper is to save computing time. The method lets you take a bigger timesteps when the curve is smoother. E.g. you take half of the time step (delta k), to evaluate the time. If you try to divide a function in different time steps, and try to approximate te slope. When it's very flat, you take double the steps, to save time and computing.
Here's a Runge-Kutta explainer video and a good resource on getting started with Python, the book Python Programming And Numerical Methods: A Guide For Engineers And Scientists.
For solving this equation, I'm assuming: $$z(t)=A \cos \omega t$$ With: $$z(0) = A$$ $$\dot{z}(0)=0$$ $$\theta(0)=0 $$ $$\dot{\theta}=0$$
Where: $A$ = amplitude $\theta$ = angular displacement $\dot\theta$ = angular velocity
Steps to work through this problem:
Python ODE solvers: https://pythonnumericalmethods.berkeley.edu/notebooks/chapter22.06-Python-ODE-Solvers.html
Simulation with matplotlib: https://towardsdatascience.com/animations-with-matplotlib-d96375c5442c
Oscillation of a simple pendulum: https://www.acs.psu.edu/drussell/Demos/Pendulum/Pendulum.html
The first step is to create a basic scatterplot with Matplotlib, and play with the positions of the points to develop intuition for the code.
import matplotlib.pyplot as plt
x = [5,7,10,7,2,17,2,9,4,11,12,9,6]
y = [99,86,95,88,111,86,103,87,94,78,77,85,86]
plt.scatter(x, y)
plt.show()
Next, draw a line between two dots by adding plt.plot(x,y). This starts to look like a pendulum.
import matplotlib.pyplot as plt
import numpy as np
x=[1,1]
y=[1.5,0.5]
plt.scatter(x,y)
plt.plot(x,y)
plt.show()
Tiny side path: I'd like to learn how to represent multiple plots in one diagram. According to the matplotlib reference guide this can be done with plt.subplots and axes created as arrays. More background info on the structure and hierarchy can be found in this tutorial about artist objects. Another method is to use nested gridspecs but that's for another time. Let's try it out.
import matplotlib.pyplot as plt
import numpy as np
# Insert data of both plots
x1 = [5,7,10,7,2,17,2,9,4,11,12,9,6]
y1 = [99,86,95,88,111,86,103,87,94,78,77,85,86]
x2 = [1,1]
y2 = [1.5,0.5]
fig, (ax1, ax2) = plt.subplots(2, 1) # Indicate how many subplots there are and their position (rows, columns)
fig.suptitle('Experiment with vertically stacked subplots')
# First plot
ax1.scatter(x1, y1)
ax1.set_ylabel('Scatterplot') # Add vertical label
# Second plot
ax2.scatter(x2, y2)
ax2.plot(x2, y2)
ax2.set_ylabel('Pendulum starting point')
plt.show()
Next up: horizontal stacking. I'm adding a constrained layout to avoid vertical labels to overlap other plots.
import matplotlib.pyplot as plt
import numpy as np
# Insert data of both plots
x1 = [5,7,10,7,2,17,2,9,4,11,12,9,6]
y1 = [99,86,95,88,111,86,103,87,94,78,77,85,86]
x2 = [1,1]
y2 = [1.5,0.5]
fig, (ax1, ax2) = plt.subplots(1, 2, layout='constrained') # Adding layout to avoid the vertical label of right plot to overlap left plot.
fig.suptitle('Experiment with horizontally stacked subplots')
# First plot
ax1.scatter(x1, y1)
ax1.set_ylabel('Scatterplot') # Add vertical label
# Second plot
ax2.scatter(x2, y2)
ax2.plot(x2, y2)
ax2.set_ylabel('Pendulum starting point')
plt.show()
Using the animation function in Matplotlib we can create moving visuals. For inspiration, here is an animated example of a pendulum. We start with trying to get a simpler example of a wave function to work. This animation uses Ipython. I installed ffmpeg with the help of this blog to play the animation.
# Credits: https://jakevdp.github.io/blog/2012/08/18/matplotlib-animation-tutorial/
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython import display
fig, ax = plt.subplots()
x = np.arange(0, 2*np.pi, 0.01)
line, = ax.plot(x, np.sin(x))
def animate(i):
line.set_ydata(np.sin(x + i / 50)) # update the data.
return line,
anim = animation.FuncAnimation( # changed 'ani' into 'anim'
fig, animate, interval=20, blit=True, save_count=50)
# To save the animation, use e.g.
#
anim.save("movie.mp4")
#
# or
#
# writer = animation.FFMpegWriter(
# fps=15, metadata=dict(artist='Me'), bitrate=1800)
# ani.save("movie.mp4", writer=writer)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
plt.close()
This works! Let's move on to creating a simple pendulum. The period of a pendulum is:
$$T=2 \pi \sqrt{\frac{L}{g}}$$The simple harmonic solution is:
$$z(t)=A \cos \omega t$$Example
Hamiltonian = total energy in the system (kinetic + potential)
$L \sin \theta$ = what we care about for potential energy
Potential energy = $mgh$ = mass * gravity * height (=$sin\theta$)
$L \sin \theta m g $
kinetic energy = $\frac{1}{2}mv^2$
theta dot is rate of energy
This youtube tutorial.
Note: Lagrange & Hamiltonian mechanics as a shortcut to the equations
Math versions of what will be coded in Python.
h is tiny small step, dt/dx.
$$x_{n+1}=x_n + h y_n $$$$y_{n+1}=y_n - \frac{hg}{l} \sin x_n $$
Can be written as:
$$\theta_{n+1} = \theta _n + h \dot\theta_n$$$$\dot\theta_{n+1} = \dot\theta _n - \frac{hg}{l} \sin (\theta _n)$$pendulum drawing
$$x=l \sin \theta$$$$y= -l \cos \theta $$import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import math
from IPython import display
# Constants
mass = 1 # Mass (kg)
length = 1 # Length of pendulum (m)
gravity = 9.81 # Acceleration/gravity (m.s-2)
ω = 10 # Natural frequency of the motion
A = 1 # Amplitude
st = 0 # Starting time (sec)
et = 15 # Ending time (sec)
time_step = 0.01 # Time step (sec)
# Initial conditions of the pendulum
theta = np.pi/4
theta_dot = 0
# other way to write this
# theta0, v0 = np.radians(60), 0 # Initial angular displacement and velocity
# Euler integration
def euler(theta, theta_dot):
new_theta = theta + time_step * theta_dot
new_theta_dot = theta_dot - time_step * gravity / length * np.sin(theta)
return new_theta, new_theta_dot
# Indentation is crucial in Python!! Not like html where it doesn't matter { }
# Python knows the difference between a float and integer. Make sure to
for _ in range(200): # range is the amount of outputs (underscore is an unused variable)
theta, theta_dot = euler(theta, theta_dot)
print(theta)
0.7853981633974483 0.7847044916451043 0.7833171481404162 0.7812366142307362 0.7784638536111858 0.7750003139867274 0.7708479293836126 0.766009123090526 0.7604868112023798 0.7542844067323942 0.7474058242508673 0.7398554850019046 0.7316383224423774 0.7227597881405564 0.7132258579652365 0.703043038489794 0.6922183735295306 0.6807594507249166 0.6686744080779959 0.6559719403443268 0.6426613051784491 0.6287523289270781 0.6142554119610675 0.599181533434752 0.5835422553596208 0.5673497258784687 0.5506166816262754 0.5333564490651458 0.5155829446827467 0.4973106739468654 0.47855472891300843 0.45933078438740493 0.43965509255438007 0.41954447598482725 0.3990163189514213 0.37808855698624705 0.3567796646276208 0.3351086413149943 0.3130949954038623 0.2907587262864473 0.2681203046184843 0.24520065066753885 0.22202111081380135 0.198603432250051 0.17496973594328152 0.15114248793614424 0.1271444690816911 0.10299874331969197 0.07872862461685556 0.054357642706406764 0.029909507774476388 0.0054080742514667005 -0.019122696124192835 -0.043658771794832074 -0.06817608924386553 -0.09265059104252901 -0.1170582638956552 -0.14137517649878995 -0.16557751702055293 -0.18964163002765655 -0.21354405267536514 -0.23726154999332624 -0.2607711491055302 -0.28405017223352624 -0.30707626834379276 -0.32982744331315245 -0.35228208850016485 -0.3744190076253135 -0.3962174418783382 -0.4176570931870372 -0.43871814559807154 -0.45938128473654755 -0.4796277153272387 -0.4994391767760523 -0.5187979568255815 -0.5376869033131539 -0.5560894340735629 -0.5739895450415219 -0.5913718166207261 -0.6082214183971546 -0.6245241122838444 -0.640266254192775 -0.6554347943366977 -0.6700172762697194 -0.6840018347802255 -0.6973771927533139 -0.7101326571223605 -0.7222581140306858 -0.7337440233246016 -0.7445814124984468 -0.7547618702106424 -0.7642775394873742 -0.7731211107273215 -0.7812858146169669 -0.7887654150615079 -0.7955542022313257 -0.8016469858184057 -0.8070390885911173 -0.8117263403294014 -0.8157050722157327 -0.8189721117502774 -0.8215247782514846 -0.8233608789959841 -0.8244787060441318 -0.8248770337898917 -0.8245551172659812 -0.8235126912273609 -0.8217499700282477 -0.8192676482998678 -0.8160669024281823 -0.8121493928228061 -0.8075172669603327 -0.8021731631772696 -0.7961202151798203 -0.7893620572298133 -0.7819028299582208 -0.7737471867499398 -0.7649003006358626 -0.7553678716207749 -0.7451561343683294 -0.7342718661572909 -0.7227223950164945 -0.7105156079395462 -0.6976599590743002 -0.6841644777766247 -0.6700387764129961 -0.6552930577921159 -0.6399381221021002 -0.6239853732269366 -0.6074468243139094 -0.5903351024626562 -0.5726634524065018 -0.5544457390578068 -0.5356964487913273 -0.5164306893430797 -0.4966641882069785 -0.4764132894176187 -0.45569494861502013 -0.4345267262959515 -0.4129267791665947 -0.3909138495227671 -0.3685072525966329 -0.3457268618227371 -0.32259309199118347 -0.2991268802717365 -0.2753496651094091 -0.25128336300954207 -0.22695034324830454 -0.20237340056274117 -0.17757572589275067 -0.15258087526546932 -0.12741273693021823 -0.10209549686921399 -0.07665360282540014 -0.05111172700379858 -0.025494727616480802 0.0001723905465878979 0.025864516328162872 0.05155647299461248 0.07722305939943581 0.10283909130742976 0.12837944266542786 0.15381908660569335 0.17913313597047212 0.2042968831508663 0.2292858390400056 0.254075770909369 0.27864273902791037 0.3029631318562025 0.32701369966195 0.35077158641871903 0.37421435986637075 0.3973200396292212 0.4200671233061479 0.4424346104654686 0.46440202449619195 0.4859494322859482 0.5070574617143205 0.5277073169682112 0.5478807917030951 0.5675602800903646 0.586728785806308 0.6053699290324597 0.6234679515500124 0.641007720022615 0.6579747275721307 0.6743550937607846 0.6901355631005521 0.705303502216666 0.7198468957967548 0.7337543414604227 0.7470150436860971 0.759618806932759 0.7715560280938275 0.7828176884190534 0.7933953450378887 0.8032811222145206 0.8124677024606773 0.8209483176275137 0.8287167400924619
What we see here, is that the Euler integration, The amplitude is getting bigger in the minus area. This should not be the case. Verified that something reasonable gets out. Let's plot it.
Instead of overridig them with _ we want to store them. The next step is to
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import math
from IPython import display
# Constants
mass = 1 # Mass (kg)
length = 1 # Length of pendulum (m)
gravity = 9.81 # Acceleration/gravity (m.s-2)
ω = 10 # Natural frequency of the motion
amplitude = 1 # Amplitude
st = 0 # Starting time (sec)
et = 15 # Ending time (sec)
time_step = 0.01 # Time step (sec)
# Initial conditions of the pendulum
theta = np.pi/4
theta_dot = 0
# other way to write this
# theta0, v0 = np.radians(60), 0 # Initial angular displacement and velocity
# Euler integration
def euler(theta, theta_dot):
new_theta = theta + time_step * theta_dot
new_theta_dot = theta_dot - time_step * gravity / length * np.sin(theta)
return new_theta, new_theta_dot
# Indentation is crucial in Python!! Not like html where it doesn't matter { }
# Python knows the difference between a float and integer. Make sure to
# Make two lists
all_thetas = [theta]
all_thetas_dot = [theta_dot]
for _ in range(200): # range is the amount of outputs (underscore is an unused variable)
theta, theta_dot = euler(theta, theta_dot)
all_thetas.append(theta) # append adds an element to a list
all_thetas_dot.append(theta_dot)
# Each list should have 101 outputs.
print("length of list:", len(all_thetas)) #len shows the length of the list
length of list: 201
That's exaclty what we want! Next, we want to create the same for time, to make sure they line up. Otherwise no plot will be generated. Here is the plot:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import math
from IPython import display
# Constants
mass = 1 # Mass (kg)
length = 1 # Length of pendulum (m)
gravity = 9.81 # Acceleration/gravity (m.s-2)
time_step = 0.01 # Time step (sec)
n_steps = 1000 # Number of steps
ω = 10 # Natural frequency of the motion
amplitude = 1 # Amplitude
# Initial conditions of the pendulum
theta = np.pi/4
theta_dot = 0
# other way to write this
# theta0, v0 = np.radians(60), 0 # Initial angular displacement and velocity
# Euler integration
def euler(theta, theta_dot):
new_theta = theta + time_step * theta_dot
new_theta_dot = theta_dot - time_step * gravity / length * np.sin(theta)
return new_theta, new_theta_dot
# Indentation is crucial in Python!! Not like html where it doesn't matter { }
# Python knows the difference between a float and integer. Make sure to
# Make two lists
all_thetas = [theta]
all_thetas_dot = [theta_dot]
# Make a list comprehension for time
all_times = [i* time_step for i in range(n_steps + 1)]
for _ in range(n_steps): # range is the amount of outputs (underscore is an unused variable)
theta, theta_dot = euler(theta, theta_dot)
all_thetas.append(theta) # append adds an element to a list
all_thetas_dot.append(theta_dot)
# Each list should have 101 outputs.
print("length of list:", len(all_thetas)) #len shows the length of the list
print("length of list:", len(all_times)) #len shows the length of the list
fig, ax = plt.subplots()
# x coordinates are time, y theta
line, = ax.plot(all_times, all_thetas)
fig.suptitle('Euler method with n = 1000 steps')
ax.set_ylabel('theta')
ax.set_xlabel('time')
plt.show()
length of list: 1001 length of list: 1001
As we noticed with the numerical output above, the amplitude of the pendulum's movement increases with each oscillation. This shows why the Euler method is not so accurate for this purpose. Yet this plot still seems relatively linear. To make this even more clear, let's increase the number of steps from 1000 to 5000.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import math
from IPython import display
# Constants
mass = 1 # Mass (kg)
length = 1 # Length of pendulum (m)
gravity = 9.81 # Acceleration/gravity (m.s-2)
time_step = 0.01 # Time step (sec)
n_steps = 5000 # Number of steps
ω = 10 # Natural frequency of the motion
amplitude = 1 # Amplitude
# Initial conditions of the pendulum
theta = np.pi/4
theta_dot = 0
# other way to write this
# theta0, v0 = np.radians(60), 0 # Initial angular displacement and velocity
# Euler integration
def euler(theta, theta_dot):
new_theta = theta + time_step * theta_dot
new_theta_dot = theta_dot - time_step * gravity / length * np.sin(theta)
return new_theta, new_theta_dot
# Indentation is crucial in Python!! Not like html where it doesn't matter { }
# Python knows the difference between a float and integer. Make sure to
# Make two lists
all_thetas = [theta]
all_thetas_dot = [theta_dot]
# Make a list comprehension for time
all_times = [i* time_step for i in range(n_steps + 1)]
for _ in range(n_steps): # range is the amount of outputs (underscore is an unused variable)
theta, theta_dot = euler(theta, theta_dot)
all_thetas.append(theta) # append adds an element to a list
all_thetas_dot.append(theta_dot)
fig, ax = plt.subplots()
# x coordinates are time, y theta
line, = ax.plot(all_times, all_thetas)
fig.suptitle('Euler method with n = 5000 steps')
ax.set_ylabel('theta')
ax.set_xlabel('time')
plt.show()
What we see now, is that the pendulum completely spirals out of control. It stops oscillating and starts spinning.
Now that the basics are covered, we're adding in the z. It is given that we should assume periodic motion.
For solving this equation, I'm assuming: $$z(t)=A \sin (\omega t)$$ $$\dot z(t)=A\omega \cos (\omega t)$$ $$\ddot z(t)=-A \omega ^2 \sin (\omega t)$$
The given equation is:
$$l \ddot{\theta} + (g+\ddot{z} ) \sin \theta = 0 $$Plugging in:
$$l \ddot{\theta} + \left( g -A \omega ^2 \sin (\omega t) \right) \sin \theta = 0 $$Make this:
$$\ddot \theta (t) = - \frac{1}{l} \left( g -A \omega ^2 \sin (\omega t) \right) \sin (\theta (t)) = 0 $$Now do Euler, because it's a second derivative, we express it in two first-order equations.
Euler: $$ \theta _{n+1}=\theta_n + h \dot\theta_n$$
$$ \dot\theta _{n+1}= \dot\theta_n - \frac{h}{l} \left( g -A \omega ^2 \sin (\omega t) \right) + sin (\theta _n)$$Then write these equations in python.
It becomes a long equation. You can always make a temporary variable. E.g:
tmp = g - A * w * w * p.sin(..)
new_theta_dot = theta_dot + h / l * tmp * ...
Amplitude & frequency - show what you get over time.
With: $$z(0) = A$$ $$\dot{z}(0)=0$$ $$\theta(0)=0 $$ $$\dot{\theta}=0$$
Where: $A$ = amplitude $\theta$ = angular displacement $\dot\theta$ = angular velocity
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import math
from IPython import display
# Constants
mass = 1 # Mass (kg)
length = 1 # Length of pendulum (m)
gravity = 9.81 # Acceleration/gravity (m.s-2)
time_step = 0.01 # Time step (sec)
n_steps = 5000 # Number of steps
omega = 10 # Natural frequency of the motion ω
A = 1 # Amplitude
# Initial conditions of the pendulum
theta = np.pi/4
theta_dot = 0
tmp = gravity - A * omega * omega * np.sin(omega * time_step) # Assuming t from the equation is time step?
# Input the new equation with z here
def euler(theta, theta_dot):
new_theta = theta + time_step * theta_dot
new_theta_dot = theta_dot - time_step / length * tmp * np.sin(theta)
return new_theta, new_theta_dot
# Make two lists
all_thetas = [theta]
all_thetas_dot = [theta_dot]
# Make a list comprehension for time
all_times = [i* time_step for i in range(n_steps + 1)]
for _ in range(n_steps): # range is the amount of outputs (underscore is an unused variable)
theta, theta_dot = euler(theta, theta_dot)
all_thetas.append(theta) # append adds an element to a list
all_thetas_dot.append(theta_dot)
# Each list should have 101 outputs.
print("length of list:", len(all_thetas)) #len shows the length of the list
print("length of list:", len(all_times)) #len shows the length of the list
fig, ax = plt.subplots()
# x coordinates are time, y theta
line, = ax.plot(all_times, all_thetas)
fig.suptitle('Euler method with 5000 steps')
ax.set_ylabel('theta')
ax.set_xlabel('time')
plt.show()
length of list: 5001 length of list: 5001
Is this correct? Does t in the equation correspond with time_step?
Next, let's move to Runge Kutta for a more accurate wave.
I removed the SciPi code in the example below, because I don't have that package installed yet.
# Example from: https://scipython.com/book2/chapter-7-matplotlib/problems/p77/animating-a-pendulum/
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# Bob mass (kg), pendulum length (m), acceleration due to gravity (m.s-2).
m, L, g = 1, 1, 9.81
# Initial angular displacement (rad), tangential velocity (m.s-1)
theta0, v0 = np.radians(60), 0
# Estimate of the period using the harmonic (small displacement) approximation.
# The real period will be longer than this.
Tharm = 2 * np.pi * np.sqrt(L / g)
# Time step for numerical integration of the equation of motion (s).
dt = 0.001
# Initial angular position and tangential velocity.
theta, v = [theta0], [v0]
old_theta = theta0
i = 0
while True:
# Forward Euler method for numerical integration of the ODE.
i += 1
t = i * dt
# Update the bob position using its updated angle.
old_theta, old_v = theta[-1], v[-1]
omega = old_v / L
new_theta = old_theta - omega * dt
# Tangential acceleration.
acc = g * np.sin(old_theta)
# Update the tangential velocity.
new_v = old_v + acc * dt
if t > Tharm and new_v * old_v < 0:
# At the second turning point in velocity we're back where we started,
# i.e. we have completed one period and can finish the simulation.
break
theta.append(new_theta)
v.append(new_v)
# Calculate the estimated pendulum period, T, from our numerical integration,
# and the "exact" value in terms of the complete elliptic integral of the
# first kind.
nsteps = len(theta)
T = nsteps * dt
print('Calculated period, T = {} s'.format(T))
print('Estimated small-displacement angle period, Tharm = {} s'.format(Tharm))
k = np.sin(theta0/2)
print('SciPy calculated period, T = {} s'
.format(2 * Tharm / np.pi * (k**2)))
def get_coords(th):
"""Return the (x, y) coordinates of the bob at angle th."""
return L * np.sin(th), -L * np.cos(th)
# Initialize the animation plot. Make the aspect ratio equal so it looks right.
fig = plt.figure()
ax = fig.add_subplot(aspect='equal')
# The pendulum rod, in its initial position.
x0, y0 = get_coords(theta0)
line, = ax.plot([0, x0], [0, y0], lw=3, c='k')
# The pendulum bob: set zorder so that it is drawn over the pendulum rod.
bob_radius = 0.08
circle = ax.add_patch(plt.Circle(get_coords(theta0), bob_radius,
fc='r', zorder=3))
# Set the plot limits so that the pendulum has room to swing!
ax.set_xlim(-L*1.2, L*1.2)
ax.set_ylim(-L*1.2, L*1.2)
def animate(i):
"""Update the animation at frame i."""
x, y = get_coords(theta[i])
line.set_data([0, x], [0, y])
circle.set_center((x, y))
nframes = nsteps
interval = dt * 1000
anim = animation.FuncAnimation(fig, animate, frames=nframes, repeat=True,
interval=interval)
plt.show()
Calculated period, T = 2.1550000000000002 s Estimated small-displacement angle period, Tharm = 2.0060666807106475 s SciPy calculated period, T = 0.3192754284070504 s
Here is an example of how to create elapsed time in Python on Stackoverflow.
import time
start = time.time()
print("hello")
end = time.time()
print(end - start)
hello 0.0002880096435546875
Here is an example of a double pendulum, here another, and a third.
