Simulating a 1D Wave on a String Using Finite Differences

The goal of this notebook is to implement a simple 1D wave propagation simulation, as explained by Hans-Petter Langtangen in his class INF5620 and more precisely in this PDF. The problem to be solved is that of a wave propagating on a string, just like in the case of a guitar.

In this notebook, we will use the simulation code writte by Mr Langtangen for the case of a (transversely) vibrating string and perform extractions of the simulated date so as to visualize it. In this process, we will use Matplotlib to embed animations and perform a little bit of Fourier analysis.

First, let's simulate the string!

Simulation¶

Below, I've adapted the code of Mr Langtangen (from the above cited PDF).

In [1]:
from numpy import linspace, zeros

def solver(I, V, f, c, L, Nx, C, T, user_action=None):
"""Solve u_tt=c^2*u_xx + f on (0,L)x(0,T]."""
x = linspace(0, L, Nx+1) # Mesh points in space
dx = x[1] - x[0]
dt = C*dx/c
Nt = int(round(T/dt))
t = linspace(0, Nt*dt, Nt+1) # Mesh points in time
C2 = C**2 # Help variable in the scheme

if f is None or f == 0 :
f = lambda x, t: 0

if V is None or V == 0:
V = lambda x: 0
u = zeros(Nx+1) # Solution array at new time level
u_1 = zeros(Nx+1) # Solution at 1 time level back
u_2 = zeros(Nx+1) # Solution at 2 time levels back

import time; t0 = time.clock() # for measuring CPU time

# Load initial condition into u_1
for i in range(0,Nx+1):
u_1[i] = I(x[i])
if user_action is not None:
user_action(u_1, x, t, 0)

# Special formula for first time step
n = 0
for i in range(1, Nx):
u[i] = u_1[i] + dt*V(x[i]) + \
0.5*C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1]) + \
0.5*dt**2*f(x[i], t[n])
u[0] = 0; u[Nx] = 0

if user_action is not None:
user_action(u, x, t, 1)

# Switch variables before next step
u_2[:], u_1[:] = u_1, u

for n in range(1, Nt):
# Update all inner points at time t[n+1]
for i in range(1, Nx):
u[i] = - u_2[i] + 2*u_1[i] + \
C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1]) + \
dt**2*f(x[i], t[n])
# Insert boundary conditions
u[0] = 0; u[Nx] = 0
if user_action is not None:
if user_action(u, x, t, n+1):
break
# Switch variables before next step
u_2[:], u_1[:] = u_1, u
cpu_time = t0 - time.clock()
return u, x, t, cpu_time


The code calls a user defined routined at each simulation step.

We can use that mechanism and define a user action routine that will store the simulation snapshots, so that we can perform some later analysis.

In [2]:
def user_action(u, x, t, n):
"Defines a user action that stores the current simulated state for further usage."
global user_data
if 'x' not in user_data:
user_data['x'] = x
if 'u' not in user_data:
user_data['u'] = [(n, u)]
else:
user_data['u'].append((t[n], u.copy()))


Let's now run a simulation using initial conditions of a 1 meter long string in the shape of a quadratic function, for a duration of 1 second and a speed of sound of 100 m/s. At each time step, we keep the simulation state for later use.

In [3]:
user_data = {}
u, x, t, cpu_time = solver(lambda x: x*(1-x), None, None, 100, 1, 100, 0.1, 1., user_action=user_action)


Visualization of the data¶

Now that we have gathered our date, let's start visualizing it. We can plot the final position of the string computed by the program:

In [4]:
import matplotlib.pyplot as plt
plt.style.use('bmh')

In [5]:
%matplotlib inline

In [6]:
plt.plot(x, u)

Out[6]:
[]

Let's examine the different simulation solution steps using a widget browser:

In [7]:
from ipywidgets import interact

In [8]:
@interact
def browse_simulation_result(n=(0, len(user_data['u']) - 1)):
"Plots a frame of the simulation."
x = user_data['x']
t, u = user_data['u'][n]
plt.plot(x, u)
plt.title('t = {:.2f}'.format(t))
plt.ylim(-0.25, 0.25)
plt.xlabel('x (a. u.)')
plt.ylabel('vertical displacement (a. u.)')


We see the string oscillating and keeping its initial bell shape. On the edges of the string, the zero displacement boundary condition keeps it fixed.

Animation of the data¶

Let's now make a moviepy animation using this data. We can embed it directly in this notebook using the ipython_display function.

In [9]:
from moviepy.video.io.bindings import mplfig_to_npimage
import moviepy.editor as mpy

duration = 30

fig, ax = plt.subplots(figsize=(4, 3), facecolor='white', dpi=150)
plt.tight_layout()
def make_frame_mpl(t):
n = int(t / duration * (len(user_data['u']) - 1))
ax.clear()
x = user_data['x']
t, u = user_data['u'][n]
ax.plot(x, u)
ax.set_title('t = {:.2f}'.format(t))
ax.set_ylim(-0.25, 0.25)
ax.set_xlabel('x (a. u.)')
ax.set_ylabel('vertical displacement (a. u.)')
return mplfig_to_npimage(fig) # RGB image of the figure

animation = mpy.VideoClip(make_frame_mpl, duration=duration)
plt.close(fig)
animation.ipython_display(fps=25, width=800)

WARNING:py.warnings:C:\Anaconda3\lib\site-packages\skimage\filter\__init__.py:6: skimage_deprecation: The skimage.filter module has been renamed to skimage.filters.  This placeholder module will be removed in v0.13.
warn(skimage_deprecation('The skimage.filter module has been renamed '

100%|███████████████████████████████████████▉| 750/751 [00:51<00:00, 14.91it/s]

Out[9]: