📓 About This Tutorial

This webpage is a static view of an interactive Jupyter Notebook tutorial. To get the full interactive experience where you can run Python code, modify examples and complete exercises, click the "Open in Colab" button. This will access the Google Colab notebook from a GitHub repository github.com/dchappell2/Computational-Physics and open it in Google Colab. You will need a Google account if you want to open it in Colab.

Open In Colab

Note: You can download a pdf of the lecture slides for this tutorial here: ODE 2 slides

ODE-2: Making Integration Code Modular

Use Functions

In this section, we re-write our Euler code around two functions:

Here are the details for each function.

Derivative Function:

Euler Integration Function

🔆 Example: Exponential Growth (revised)

In this example, we re-write our exponential growth example using Euler integration and the modular approach described above. Our differential equation is $$\frac{dy}{dt}=ay.$$

Note: We want to write our code generally, so that the derivative function can take one or more parameters. In the main code we bundle the parameters in an array called params. But in the derivative function, the parameters are listed as separate variables. To handle this, we need to "unpack" the params array and assign each element of the array to each of the parameters in the derivative function. To do this we use the "unpacking" operator, i.e. we place an asterisk * before the name of our params array when we pass it to the derivative function like this:

```

f = deriv(t[n], y[n], params) # use "" to unpack params

```

Here's the code:

import numpy as np
import matplotlib.pyplot as plt

#########  Derivative Function  #########
#
# This function returns the derivative for the exponential growth ODE dy/dt = ay

#########  Derivative Function  #########

def deriv_exp(t, y, a):

    dydt = a*y
    return dydt


#########  Euler Integration  #########

def Euler(deriv, y0, tmax, dt, params):

    #########  Create Arrays  #########

    N = int(tmax/dt)+1    # number of steps in simulation
    y = np.zeros(N)       # array to store y values
    t = np.zeros(N)       # array to store times

    y[0] = y0             # assign initial value

    #########  Loop to implement the Euler update rule  #########

    for n in range(N-1):
        f = deriv(t[n], y[n], *params)  # use "*" to unpack params
        y[n+1] = y[n] + f*dt
        t[n+1] = t[n] + dt

    return t, y


#########  Parameters  #########

a    = 0.2       # decay constant
tmax = 100       # maximum time
dt   = 0.5       # time step
y0   = 1         # initial value of y

params = [a]     # bundle parameters in an array


#########  Perform Euler Integration  #########

t, y =  Euler(deriv_exp, y0, tmax, dt, params)


#########  Analytic Solution  #########

y_true = y0 * np.exp(a*t)


#########  Plot Solution  #########

plt.plot(t, y, label='Euler')
plt.plot(t, y_true, '--', label='Analytic')

plt.xlabel('t')                  # label the x and y axes
plt.ylabel('y')
plt.title("Exponential Decay")   # give the plot a title
plt.legend()                     # display the legend
plt.show()                       # display the plot
Output image

✅ Skill Check 1

In this example, you will solve for the number of photons $n$ in a laser cavity. Use the example above to solve the following ODE for $n(t)$: $$ \frac{dn}{dt} = (GN_0 -k) n- (\alpha G) n^2$$ where $n$ is the number of photons in the laser field and

Specifications:

Run your code to produce two plots:

###  Your Code Here

Use Arrays to Handle Systems of ODEs

We've already seen one system of ODEs in the previous chapter: the simple hamonic oscillator.

Second-order ODE's arise naturally from Newton's second law since acceleration is a second order derivative:

$$\frac{d^2x}{dt^2}=a = \frac{F}{m}$$

where the acceleration is potentially a function of position, velocity, time and mass: $a=F(x, v, t)/m$.

As discussed in the lecture slides, we write this second-order equation as two first-order equations like this:

$$ \frac{dx}{dt} = v $$

$$ \frac{dv}{dt} = \frac{a}{m}. $$

We can write this system of ODEs as a single vector equation, where the components of the generalized vector are position and velocity. If we define our generalized vector as $\vec{y}$, its components will be $$y_0(t) = x(t)$$ and $$y_1(t) = v(t).$$ The vector $\vec{y}$ may be written as

$$\vec{y} = \begin{pmatrix}y_0(t) \\ y_1(t) \end{pmatrix} = \begin{pmatrix}x \\ v \end{pmatrix}.$$

With these substitutions, we can write our system of coupled ODEs as

$$ \frac{dy_0}{dt} = y_1 $$

$$ \frac{dy_1}{dt} = -(k/m)y_0. $$

Our system of differential equations may now be written as a single vector equation:

$$\frac{d\vec{y}}{dt}=\vec{a}(\vec{y},t),$$ where $\vec{a}(\vec{y},t)$ is a vector-valued acceleration:

$$\vec{a}(\vec{y},t)= \begin{pmatrix}v \\ -(k/m)x \end{pmatrix} = \begin{pmatrix}y_1 \\ -(k/m)y_0\end{pmatrix}.$$

We can now apply our favoriate numerical integration method to this single vector-valued ODE. For example the Euler update rule for our generalized state vector would be:

$$ \vec{y}_{n+1} = \vec{y}_n + \vec{a}_n \Delta t.$$

🔆 Example: Vectorized SHO example

Here's the code for the Simple Harmonic Oscillator. We introduce a new derivative function and a new Euler function that handles multi variables.

import numpy as np
import matplotlib.pyplot as plt

#########  Derivative Function  #########
#
# This function returns the derivatives for the
# Simple Harmonic Oscillator ODE

def deriv_sho(t, y, m, k):

    # extract variables from y array
    x = y[0]          # position
    v = y[1]          # velocity

    # calculate derivatives
    dxdt = v
    dvdt = -k/m*x

    # return derivatives in a numpy array
    return np.array([dxdt, dvdt])


#########  Euler Integration  #########
#
# This function integrates an ODE specified by a passed derivative function
#
# Passed parameters:
#   deriv  = function created by the user that returns the derivatives for the system of ODEs
#   y0     = 1xN array containing initial conditions
#   tmax   = maximum time to run the integration
#   dt     = time step
#   params = 1xM array containing parameters for the derivative function
#
# Returned values:
#   t      = 1D array containing time values of the solution
#   y      = 2D array whose columns are the solutions for each variable

#########  Multi-Variable Euler Integration  #########

def Euler_Vec(deriv, y0, tmax, dt, params):

    #########  Create Arrays  #########

    # determine the number of variables in the system from initial conditions
    nvar = 1 if not isinstance(y0, np.ndarray) else y0.size

    N = int(tmax/dt)+1         # number of steps in simulation
    y = np.zeros((N,nvar))     # array to store y values
    t = np.zeros(N)            # array to store times

    if nvar == 1:
        y[0] = y0              # assign initial value if single variable
    else:
        y[0,:] = y0            # assign vector initial values if multivariable

    #########  Loop to implement the Euler update rule  #########

    for n in range(N-1):
        f = deriv(t[n], y[n], *params)
        y[n+1] = y[n] + f*dt
        t[n+1] = t[n] + dt

    return t, y


import numpy as np
import matplotlib.pyplot as plt

#########  Parameters  #########

m    = 1         # mass
k    = 1         # spring constant
tmax = 10        # maximum time
dt   = 0.001     # time step
x0   = 1         # initial position
v0   = 0         # initial velocity

params = np.array([m,k])    # bundle parameters together
y0 = np.array([x0,v0])      # bundle initial conditions


#########  Perform Euler Integration  #########

t, y =  Euler_Vec(deriv_sho, y0, tmax, dt, params)

x = y[:,0]       # extract positions
v = y[:,1]       # extract velocities


#########  Plot Solution  #########

plt.plot(t, x, label='x')        # plot position
plt.plot(t, v, '--', label='v')  # plot velocity

plt.xlabel('t')                  # label the x and y axes
plt.ylabel('x, v')
plt.title("SHO")                 # give the plot a title
plt.legend()                     # display the legend
plt.show()                       # display the plot
Output image

✅ Skill Check 2

Write code to model the motion of a mass attached to a pivot by a spring with spring constant $k$ immersed in a fluid with fluid drag coefficient $b$. Assume the mass can move in the x-y plane.

$$\frac{dx}{dt} = v_x$$

$$\frac{dy}{dt} = v_x$$

$$\frac{dv_x}{dt} = -\frac{kx+bv_x}{m} $$

$$\frac{dv_y}{dt} = -\frac{ky+bv_y}{m} $$

### Your Code Here