There are many ways you can numerically solve differential equations, from relying on your own implementation of different methods like Euler-Cromer to using built in libraries that implement methods for you. On popular and powerful tool toward this end is Scipy's odeint tool, which I'll walk you through using in this guide. To get started we'll import the needed function as well as numpy.
import numpy as np
from scipy.integrate import odeint
Say you have any physical situation where you have found your equations of motion. Thus you have 1 or more usually 2nd order differential equations. For our example here we will take the case of a damped harmonic pendulum, but it could be any 2nd order differential equation. The equation for a damped oscillator is: $$ \ddot \phi = -\beta \dot\phi -\frac{g}{R}\sin(\phi)$$ Before worrying about the differential equation, I'm going to go ahead and define my constants here. Remember, you are solving this numerically, so you'll need numerical values to plug in for all your constants!
β = 0.2 # 1/s
g = 9.8 # m/s^2
R = 2 # m
For some convenience here, I'm also going to go ahead and define the classic angular frequency: $$ \omega = \sqrt{\frac{g}{R}} $$
ω = np.sqrt(g/R)
We know come across a very important point. Almost any numerical solver will require you to convert your 2nd order differential equation into a system of 1st order differential equations. Solvers like 1st order equations. Since most of the equations of motion you will be finding this semester will be 2nd order, you should get real familiar with the following task.
The idea is that we break the 2nd order equation up into two 1st order equations by introducing a new variable. For instance, for the above equation let's introduce a variable a, such that
$$ \dot \phi = a $$
Easy so far correct? Now, since $\dot \phi = a$, we can also say that
$$ \ddot \phi = \dot a $$
But we already know what $\ddot \phi$ is from our initial equation of motion, so we can just set $\dot a$ equal to that. Note that we will also plug in $a$ for the $\dot \phi$ that is multiplied by the $\beta$:
$$\ddot \phi = \dot a = -\beta a - \omega^2\sin(\phi)$$
Now, if you look carefully, we now have TWO first order differential equations:
\begin{align}
\dot \phi &= a \\
\dot a &= -\beta a -\omega^2\sin(\phi)
\end{align}
Important: Note that both equations have all $\dot \phi$ or $\dot a$ on the left hand side of the equation, and only $\phi$ or $a$ (or constants) on the right hand side of the equation. This is the form you'll want to make sure your equations are in before you can utilize them in a solution.
We enter this system of equations into Python in a particular format that odeint will be able to understand, by defining a function:
def DE( state, time ):
# Variables in the DE
ϕ, a = state
# RHS of equations
derivs = [a, -β*a - ω**2*np.sin(ϕ)]
return derivs
I'll break down the above a bit as it can be confusing. The function is dependent on two inputs: the current state of the system (aka, what $\phi$ and $a$ are currently) and the current time the system is at. For convenience, we first "unpack" the state to get the individual variables. These should always correspond to the variables you have on the left hand side of your system of first order differential equations (just not dotted). Then we basically define a list of the right hand side of those equations. The first entry in the list corresponds to the first variable we listed above, the second to the second listed variable, etc. So since $\dot \phi = a$, we enter in $a$ to the first spot on the list, and then enter in all of $-\beta a - \omega^2\sin(\phi)$ into the second spot. Then we return the derivatives as the output of our function. That's it!
It is important to note that your system of DE's can be however large you want here, so long as they are all first order differential equations. If you say had two equations of motion in polar coordinates, one with $\ddot r$ and one with $\ddot \phi$, you'd just break them both down into two 1st order DE's and then your system above would have a total of 4 differential equations!
We need one more thing before we can plug everything into the numerical solver. Since we know that integrating yields unknown constants, you always need to know something about the initial state of the system to be able to get a unique solution. As such, we need to dictate what our initial $\phi$ and $a$ were for the system. We'll let the pendulum start from $90^\circ$ to the vertical, and we'll let $a$ start at 0. Remember that $a=\dot \phi$, so this is the same thing as saying that we drop the pendulum from rest.
ϕi = np.pi/2
ai = 0
We also need to determine what time period we are interested in obtaining the solution over. A particular thing about solving DE's numerically is that you will only get a solution over the time interval you entered in. If you decide you want a different time interval, you need to rerun the whole solution. For this let's say we are interested in looking at the solution over the first 20 seconds of its motion. We'll create a vector containing a series of time steps from 0 to 20. We have some options in how we want to do this:
# This options lets us imput the starting and ending times as well as how many data points we want in between.
#So in this case I'm going from 0 to 20 in whatever even intervals will give me 500 data points
ts = np.linspace(0,20,500)
# This options lets us input the starting and ending times and then the step size of each time interval.
#So here I'm going from 0 to 20 by increments of 0.01 seconds.
ts = np.arange(0,20,0.01)
You only need ONE of those definitions! I just left them both in as an example. Use whichever makes the most sense to you. Now at this point we are set to obtain our solutions! odeint takes options of the form
odeint( [system of 1st order odes], [initial conditions], times, optional arguments for f)We've already defined everything we need so:
sols = odeint( DE, [ϕi, ai], ts )
If we look at the solution, we'll see that it has two columns, and each row corresponds to the state of the system at a different time. The first column will correspond to whatever your first variable in your system of DE's was (so $\phi$ for us) and the second to the second variable ($a$). At this point, you have your numeric solution! Then you'll usually want to plot it in some way to interpret your results.
I could list all the $\phi$ values:
sols[:,0]
Or list all the $a$ (aka $\dot \phi$) values:
sols[:,1]
We can check how we did by plotting things vs time. We are probably most interested in how the position of the pendulum ($\phi$) changed over time. Using your preferred plotting method:
from bokeh.io import output_notebook
from bokeh.resources import INLINE
from bokeh.plotting import figure, show
output_notebook(INLINE)
f = figure(width=800, height=400, x_axis_label='Time', y_axis_label='ϕ')
f.line(ts, sols[:,0], line_width=3, color='red')
show(f)
Looks pretty good and reasonable!