# 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]: