📓 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 1 slides

ODE-1: Euler's Method for Solving Ordinary Differential Equations (ODEs)

Euler's Method for First-Order ODE's

We developed the ideas behind Euler's method in the lecuture slides. As an overview, suppose we want to solve the following first-order ODE:

$$ \frac{dy}{dt} = f(y), $$

where $f(y)$ is some function of our variable $y$.

The Euler update rule for $y$ is a result of replacing the derivative with its discrete approximation $\frac{\Delta y}{\Delta t} \approx f(y)$. The Euler update rule is given by

$$ y_{n+1} = y_n + f(y) \Delta t.$$

The update rule for time is simply:

$$ t_{n+1} = t_n +t \Delta t.$$

🔆 Example: Constant Velocity Motion

This is a very simpl example. The equation for constant velocity motion is given by

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

where $v=$ the velocity. This formula is the simplest possible ODE. The solution is equation can be easily integrated with solution

$$x_{exact}(t) = x_0+vt.$$

Here's some Python code to implement the Euler rule:

import numpy as np
import matplotlib.pyplot as plt

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

v    = 10             # velocity
tmax = 10             # maximum time
dt   = 2              # time step
x0   = 10             # initial value of x


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

N = int(tmax/dt)+1    # number of steps in simulation

x = np.zeros(N)       # array to store positions
t = np.zeros(N)       # array to store times

x[0] = x0             # assign initial value

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

for n in range(N-1):
    x[n+1] = x[n] + v*dt   # Euler update rule for position
    t[n+1] = t[n] + dt     # update time


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

x_true = x0 + v*t


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

plt.plot(t, x, 'ro', label='Euler')
plt.plot(t, x_true, 'b--', label='Analytic')

plt.xlabel('t (s)')                   # label the x and y axes
plt.ylabel('x (m)')
plt.title("Constant Velocity Motion") # give the plot a title
plt.legend()                          # display the legend
plt.show()                            # display the plot
Output image

🔆 Example: Exponential Growth

The equation for exponential growth or decay is given by

$$ \frac{dy}{dt} = ay,$$

where $a=$ the growth (or decay) rate. This formula is of the form of our first-order ODE above wehre $f(y) = ay$. This equation can be easily integrated with solution

$$y(t) = y_0 e^{at}.$$

If $a>0$ the quantity $y$ blows up exponentially and if $a<0$ it decays to zero.

Here's some Python code to implement the Euler rule:

import numpy as np
import matplotlib.pyplot as plt

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

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


#########  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):
    y[n+1] = y[n] + a*y[n]*dt
    t[n+1] = t[n] + dt


#########  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 Growth")  # give the plot a title
plt.legend()                     # display the legend
plt.show()                       # display the plot
<Figure size 640x480 with 1 Axes>
Output image

✅ Skill Check 1

Plot the percent error produced by Euler's method using the formula

$$percent\; error=\frac{|y-y_{analytic}|}{y_{analytic}}\times 100\%$$

Label your axes and give the plot a title

# your code here

✅ Skill Check 2

Imagine numerically integrating an exponential decay problem where the growth rate $a<0$.

your answers here 

# your code here

✅ Skill Check 3

import numpy as np
import matplotlib.pyplot as plt

# your code here

Euler's Method for Second-Order ODE's

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

$$F=ma$$.

Solve for acceleration

$$a=\frac{F}{m}.$$

Write as a position derivative and explicitly show the possible position, velocity and time dependence of the force:

$$\frac{d^2x}{dt^2}=\frac{a}{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 now apply the Euler method to each of these first-order ODEs:

$$ x_{n+1} = x_n + v_n \Delta t,$$

$$ v_{n+1} = v_n + a(x_n, v_n, t)\Delta t,$$

The update rule for time:

$$ t_{n+1} = t_n +t \Delta t.$$

✅ Skill Check 4

The equation of motion for a simple harmonic oscillator is

$$\frac{d^2x}{dt^2} =-\omega_0^2 x$$

where $\omega_0^2=k/m$ is the angular frequency of the oscillations.

$$ x_{theory}(t) = x_0 \cos(\omega_0 t)$$ assuming $v_0=0$. Note: we plot the absolute error instead of the relative error because this solution passes through zero and will cause the relative error to "blow up."

# your code here

Euler-Cromer-Aspel Method

If we modify our Euler method in one shockingly simple way, we go from a first-order Euler method that is inherently unstable for the simple harmonic oscillator and other oscillatory, energy-conserving systems, to one that is stable with much smaller errors.

The simple trick is this: use the new velocity v[n+1] in the position update rule instead of the old velocity v[n]. With this one simple change, we get a much more robust integration method for oscillatory, energy-conserving systems.

The Euler-Cromer rule looks like this for the simple harmonic oscillaror:

```

for n in range(N-1):

v[n+1] = v[n] - k/mx[n]dt

x[n+1] = x[n] + v[n+1]*dt # modified to use v[n+1] instead of v[n]

t[n+1] = t[n] + dt

```

Historical Note: While this method is commonly called the Euler-Cromer or the Symplectic Euler method, we believe the addition of "Aspel" is warranted. Alan Cromer attributes the discovery of this method to a high school stutent, Abby Aspel. Richard Feynman notes that the method has been discovered, forgotten and rediscovered many times since Newton. It was Cromer's paper analyzing Aspel's re-discovery of the method, however, that came at a time (1980) when numerical computations were more widespread in physics. Cromer received recognition while the (female) discoverer did not.

✅ Skill Check 5

In this problem, we will modify our Euler method to use the Euler-Cromer-Aspel method: replace the v[n] term in the position update rule with v[n+1]. This will require movint the x update rule after the v update. The for loop with the update rules will look something like this:

```

for n in range(N-1):

v[n+1] = v[n] - k/mx[n]dt

x[n+1] = x[n] + v[n+1]*dt # modified to use v[n+1] instead of v[n]

t[n+1] = t[n] + dt

```

import numpy as np
import matplotlib.pyplot as plt

# your code here