Frolian's blog standing on the shoulders of PyGiants

• About
• Archive

Simulating a 1D Wave on a String Using Finite Differences

Thu, 21 Apr 2016

Numerical physics

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).

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.

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.

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:

import matplotlib.pyplot as plt
plt.style.use('bmh')

%matplotlib inline

plt.plot(x, u)

[<matplotlib.lines.Line2D at 0xdc4c208>]

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

from ipywidgets import interact

@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.

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]