Prelude to chaos: solving non-linear differential equations with solve_ivp and scipy
Welcome to this blog where we will use solve_ivp and Python to solve a non-linear differential equation. Non-linear differential equations often describe chaotic behaviour. In this blog we will solve a non-linear differential equation which is very interesting but not quite chaotic yet. This will help us to learn solve_ivp and scipy better and prepare for some full chaotic behaviour in a next blog.
Recap of simple harmonic motion
In my previous blog we solved the differential equation for simple harmonic oscillations using scipy, solve_ivp and Python. Remember that the equation of motion looked like
One of the systems this equation describes is the motion of a mass 𝑚 on a spring with a spring constant 𝑘. You can see that the acceleration 𝑑²𝑥/𝑑𝑡² is always in the direction of the equilibrium, 𝑥 = 0. This is a so-called restoring force. Because the mass will overshoot the equilibrium, it will start oscillating. And with no friction, it will oscillate indefinitely.
The equation for harmonic motion we covered in the previous blog can be solved analytically and the solutions are sine and cosine functions (try and plug those into the equation of motion). So, solving such an equation numerically was just an exercise. In this blog we will solve a more difficult equation that cannot be solved analytically (in other words, the solution cannot be formulated in elementary functions like sine, cosine, logarithms, etc).
Forces that make a pendulum swing
The equation we will solve today belongs to a pendulum that can rotate around a point of rotation. The pendulum we will study has all its mass at one point and is connected to the point of rotation through a rigid, mass-less, rod. In Fig. 1 you see a schematic depiction of the situation. We defined 𝜃 as the angle between the rod and the vertical at the point of rotation. We will take 𝜃 to be positive to the left.
To write down the equation of motion for the pendulum we will first need to analyse the forces acting on the pendulum. The only force that will contribute to the movement is the gravitational force (is it a force? Check my blog on curved spaces and general relativity if you are interested). The gravitational force of gravity is pointing downwards, but the component of the force in the 𝜃 direction (see Fig. 1) is what we are interested in:
You can see that it points in the negative 𝜃 direction when 𝜃 is positive and it will therefor function as a restoring force, just as with the harmonic oscillations of a mass on a spring. Thus, the solutions we find will possibly resemble harmonic motion. But let’s not get ahead of ourselves :)
The non-linear equation of motion
Before we can write down our equation of motion, we need to talk about a trick to make things more easy for ourselves. Even though our pendulum moves in two dimensions, the problem really is one dimensionally. The only free variable is the angle 𝜃. If we manage to write our equation of motion in terms of this variable, we will be of to a good start.
First we will look at the acceleration and how it translates into changes in 𝜃. The mass and the rod that holds it can only move in a circle around the point where the rod is fixed. We take 𝑥 to be the length travelled on this circle and 𝑎 the acceleration with which the mass m travels on this circle. Now we can write the following: 𝑥=𝑙 𝜃, where 𝑙 is the length of the rod. This makes sence since if we move one whole rotation we have 𝜃 = 2𝜋 and then x = 2𝜋𝑙, the circumference of a circle with radius 𝑙. We can thus write:
The equation of motion now becomes:
If we rearrange things a bit and define
we will end up with:
This equation looks very similar to the equation of motion of a mass m on a spring, with only one exception: on the right hand side we now have sin(𝜃) instead of just 𝜃. And this makes a difference!
An approximation that makes our equation of motion linear
Now we have a nice differential equation in terms of the angle 𝜃 and the time 𝑡. That we have a sin(𝜃) on the right instead of a 𝜃 makes for a non-linear differential equation. The equation of harmonic motion is linear in because the second derivative 𝑑²𝜃/𝑑𝑡² is directly proportional to the function 𝜃(t) itself. Now for our pendulum the second derivative 𝑑²𝜃/𝑑𝑡² is not proportional to 𝜃(𝑡), but to the sin(𝜃). Analytically solving non-linear differential equations is another game and more often than not, impossible analytically. But we are armed with numerical options and we will not be deterred!
Before we tackle the non-linear equation of motion we can first study it in the special case where it is linear. If we keep to small values of 𝜃, we can approximate 𝜃 ≈ sin(𝜃), which gives:
From this we learn that at small deviations we just have harmonic motion and we know that the solution of this is:
Where 𝐴 and 𝐵 are constants that are determined by the starting conditions and our solution becomes:
Lets plot some examples of this approximation using the following code:
import numpy as np
import matplotlib.pyplot as plt
# We make a function to calculate the linear version
def theta_linear(t = np.arange(0, 1, 0.001), omega = 9.8/1, th_0 = 0.01):
return t, np.array([th_0*np.cos(omega*t_i) for t_i in t])
# We plot the linear solution for three pendulum lengths
for omega in [9.8/0.5, 9.8/1.0, 9.8/2.0]:
plt.plot(*theta_linear(omega=omega), label = "$\\omega$ = "+str(omega)+" 1/s")
plt.ylim((-0.011, 0.016))
plt.legend()
plt.ylabel('$\\theta (t)$ (radians)')
plt.xlabel('time (seconds)')
In Fig. 2 you can see you get harmonics oscillations with increasing period when you increase the length of the pendulum. Lets make an interesting animation of our solution using the following code:
import matplotlib.animation as animation
from matplotlib import rc
rc('animation', html='jshtml', embed_limit=500)
def plotPendulum(t, th, l, start_at_frame = 0, interval = 100, \
ani_name = "Pendulum_small_displacement.gif"):
# We use this function to go from theta to cartesian coordinates
def cartCoords(th, l):
return -l*np.sin(th), -l*np.cos(th)
frames = len(t)-start_at_frame
# We make a figure with two subplots
fig = plt.figure(figsize=(10,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
# We add the first data points to the (x,y) plot and
# set the limits
ax1.axvline(0, color='r', linestyle='-')
x0,y0 = cartCoords(th[0], l)
ax1.set_xlim(-l-0.2, l+0.2)
ax1.set_ylim(-l-0.2, l+0.2)
ax1.set_xlabel("x (m)")
ax1.set_ylabel("y (m)")
# We add a dot and line to the (x,y) plot
dot, = ax1.plot(x0, y0, "ko")
line, = ax1.plot([0, x0], [0, y0], "k")
# We make the second plot showing theta as a
# function of time
ax2.plot(t, th, "gray")
ax2.axhline(y=0, color='r', linestyle='-')
# We also add a dot for the animation
dot2, = ax2.plot(t[0], th[0], "ko")
ax2.set_xlabel("t (sec)")
ax2.set_ylabel("$\\theta$ (radians)")
ax2.yaxis.set_label_position("right")
ax2.yaxis.tick_right()
# We can alter the start frame if we want
t = t[start_at_frame:]
th = th[start_at_frame:]
# This function updates the plot for the animation
def update(frame):
# Updating the plots
xi, yi = cartCoords(th[frame], l)
dot.set_data([xi], [yi])
line.set_data([0, xi], [0, yi])
dot2.set_data([t[frame]], [th[frame]])
# Create the animation
ani = animation.FuncAnimation(fig=fig, func=update, \
frames=frames, interval=interval)
ani.save(ani_name)
plt.close()
return ani
plotPendulum(*theta_linear(t = np.arange(0, 2.0, 0.01), omega=9.81/1, \
th_0 = 0.5), l=1, ani_name='fig3_pendulum_lin.gif')
As you can see in Fig. 3, we have nice harmonic oscillations.
Solving the non-linear equation of motion using solve_ivp and scipy
Now it is time to solve the real deal without any approximations. For readability here is the equation of motion once more:
We will go and solve this equation using scipy’s solve_ivp function. As you have learned from my previous blog solve_ivp only solves a system of first order differential equations. We need to write our second order differential equation as two coupled first order ones. Like we did in the last blog we can write the rate of change of 𝜃, lets call it the 𝜃 velocity, as the first kind of trivial equation :)
Here we used the dot one the 𝜃 to indicate the rate of change in time. Now we write the rate of change of the rate of change ;) of 𝜃 as the second equation:
We now need to code a step function that solve_ivp will call at every integration step to update the values of 𝜃 and d𝜃/dt. This function must have as first parameter the time t and as a second parameter a list containing the new values of 𝜃 and d𝜃/dt. We will call this list U and we will take the first item in the list to be 𝜃 and the second to be d𝜃/dt. This order is determined by the starting values we will give to solve_ivp. Our step function will look like this then:
def step(t, U, omega):
th, thDot = U
U_new = [thDot, -1 * omega**2 * np.sin(th)]
return U_newThe step function will return updated values of the values in U. So U[0] is 𝜃 and the step function must give the value of the time derivative of 𝜃 as an update for the integrator. So, U_new[0] will contain thDot, according to Eq. 1. And U[1] will contain d𝜃/dt, so U_new[1] must contain the time derivative of d𝜃/dt, which is −𝜔²sin(𝜃). Summarizing all this gives:
So the most difficult part of coding solve_ivp is getting the eq. 1 and eq. 2 coded in the 𝑈 and 𝑈_new like above. If you want to see some more examples, check out my previous blog. Next, we will make a function that implements solve_ivp and solves the non-linear differential equation:
from scipy.integrate import solve_ivp
def theta_non_linear(th0, thDot0, omega=9.81/1, \
curve_start = 0, curve_end = 3, nr_points_curve = 100):
# Starting conditions for theta (th0) and the time derivative of
# theta (thDot0)
U0 = [th0, thDot0]
# The time values for which we want out solution
t_pts = np.linspace(curve_start, curve_end, nr_points_curve)
# Start the integrator
result = solve_ivp(step, (curve_start, curve_end), U0, \
t_eval=t_pts, args=(omega,))
# Return the time values and corresponding theta values (result.y has
# the solutions for theta and thetaDot in the same order as you
# have given them in the starting conditions)
return result.t, result.y[0]Exploring the solutions of solve_ivp in the non-linear case
Let’s have a look at the different possible solutions we can find for the non-linear case and how it compares to the approximate solutions of the linear case. The function below plots four cases for four different starting conditions. With different starting conditions I mean different values of 𝜃_0, so different starting displacements. I keep the starting velocity of the pendulum the same, namely at zero (you just let loose the pendulum at an angle 𝜃_0).
def make_compare_plot():
fig = plt.figure(figsize=(10,10))
ax1 = fig.add_subplot(221)
ax2 = fig.add_subplot(222)
ax3 = fig.add_subplot(223)
ax4 = fig.add_subplot(224)
ax1.set_title('$\\theta_0 = 0.05 \pi$')
ax1.plot(*theta_linear(t = np.arange(0, 3.0, 0.01), omega=9.81/1, th_0 = 0.05), "k", label="Linear")
ax1.plot(*theta_non_linear(th0 = 0.05, thDot0=0.0, omega=9.81/1), "r.", label="Non-linear")
ax1.set_ylabel("y (m)")
ax2.set_title('$\\theta_0 = 1/2 \pi$')
ax2.plot(*theta_linear(t = np.arange(0, 3.0, 0.01), omega=9.8/1, th_0 = 1/2*np.pi), "k")
ax2.plot(*theta_non_linear(th0 = 1/2*np.pi, thDot0=0.0, omega=9.8/1), "r.")
ax2.set_ylabel("$\\theta$ (radians)")
ax2.yaxis.set_label_position("right")
ax2.yaxis.tick_right()
ax3.set_title('$\\theta_0 = 0.9 \pi$')
ax3.plot(*theta_linear(t = np.arange(0, 3.0, 0.01), omega=9.8/1, th_0 = 0.9*np.pi), "k")
ax3.plot(*theta_non_linear(th0 = 0.9*np.pi, thDot0=0.0, omega=9.8/1, curve_end=3, nr_points_curve = 100), "r.")
ax3.set_xlabel("x (m)")
ax3.set_ylabel("y (m)")
ax4.set_title('$\\theta_0 = \pi$')
ax4.plot(*theta_linear(t = np.arange(0, 30.0, 0.01), omega=9.8/1, th_0 = np.pi), "k")
ax4.plot(*theta_non_linear(th0 = np.pi, thDot0=0.0, omega=9.8/1, curve_end=30, nr_points_curve = 1000), "r.")
ax4.set_xlabel("t (sec)")
ax4.set_ylabel("$\\theta$ (radians)")
ax4.yaxis.set_label_position("right")
ax4.yaxis.tick_right()
In Fig. 4 you see that at 𝜃_0 = 0.05𝜋 you see there is no difference between the linear and non-linear case. This is understandable because the approximation 𝜃 ≈ sin(𝜃) that we used to get a linear differential equation, is valid at small 𝜃. When we increase 𝜃_0 to 1/2𝜋 we see our approximation breaking down and the linear equation is no longer a good description of the system. At 𝜃_0 = 0.9𝜋 we see the difference even stronger.
Now the last case in Fig. 4 we put the pendulum upside down at a displacement of 𝜃_0 = 𝜋 . This is a funny case. In this position, it should stay there since there is no force pushing it to the left or right. But it is not a stable solution, only a minute small deviation from exactly 𝜃_0 = 𝜋 and the pendulum will fall and start rotating. So a small rounding error or another numerical effect can throw it off.
If we run:
plotPendulum(*theta_non_linear(th0=np.pi, thDot0=0, omega=9.8/1, \
curve_start = 0, curve_end = 20, nr_points_curve = 500), \
l = 1, start_at_frame = 60, interval = 20)we get the animation of Fig. 5. In this figure you can see an animation of this funny case of an upside down pendulum. For our calculations it is just a minute little numerical deviation from 𝜃_0 = 𝜋 that will throw it off and the rotation starts.
An exercise for you
As an exercise for you to get a better feel for the non-linear behaviour of this system, solve the equation for (𝜃_0 = 𝜋) but use very small deviations to the length of the pendulum. For example, for l = 1.0, l = 1.01, l=1.001 and l=1.002. What do you see? Do such small deviations make any impact on the solution? Let me know in the comments below what you find!
Also, are you working on some fun project for which you need solve_ivp? Let me know in the comments, I would love to know what you are brewing!
If you’ve enjoyed this post, clap 👏 for the story, leave a comment and follow me for more exciting physics and tutorials on Python and numerical techniques.
Take care and cheers,
Ben
ps: In a next blog we will use solve_ivp to calculate trajectories of photons around black holes. In preparation you can already have a look at one of my previous blogs on Calculating lengths in curved spaces using SymPy’s symbolic mathematics, Python and Matplotlib.
