Solving the tensioned string equation in time domain and frequency domain

Waves Physics Numerical Physics Animations MOOC

In this blog post, we will look at the vibrating string equation and solve it using two different mathematical methods. First, we will use a time-domain finite difference method. In a second step, we will use a frequency domain approach (i.e. a Fourier method) to solve the same problem.

I wrote this post while going through the excellent MOOC Fundamentals of Waves and Vibrations. Feel free to check out what I wrote about the experience here.

The vibrating string equation

Let's start by introducing the vibrating string equation.

The equation of the transverse displacement of a vibrating string $y(x,t)$ is, $\forall x \in [0, L]$, $\forall t \geq 0$:

$$ m \dfrac{\partial^2 y}{\partial t^2} - T \dfrac{\partial^2 y}{\partial x^2} = 0 $$

We will solve this with so-called Dirichlet boundary conditions, meaning that the string has a zero displacement at its ends, i.e. $\forall t \geq 0 \; y(0, t) = y(L, t) = 0$. Initial conditions will be as follows:

  • $\forall x \in [0, L], \quad y(x, 0) = y_0(x)$
  • $\forall x \in [0, L], \quad \dot{y}(x, 0) = 0$

Time-domain solution

Let's implement a time-domain solution known as finite differences. The partial differential equation above, once discretized on a 1D grid, becomes:

$$ m \dfrac{y_i^{n+1} - 2 y_i^{n} + y_i^{n-1}}{\Delta t^2} - T \dfrac{y_{i+1}^{n} - 2 y_{i}^{n} + y_{i-1}^{n}}{\Delta x^2} = 0 $$

After reordering the terms, we obtain the following propagation step:

$$ \forall i \in [1, N-1] \quad y_i^{n+1} = 2 y_i^{n} - y_i^{n-1} + \dfrac{T}{m} \dfrac{\Delta t^2}{\Delta x^2} \left ( y_{i+1}^{n} - 2 y_{i}^{n} + y_{i-1}^{n}\right ) $$

We can setup a simple implementation of the solution as follows:

In [1]:
import numpy as np

class TimeDomainSolver:
    
    def __init__(self, L, nx, T, m, y0):
        """Setup of the time domain solver."""
        self.x = np.linspace(0, L, num=nx)
        self.dx = self.x[1] - self.x[0]
        c = np.sqrt(T/m)
        self.dt = self.dx / c
        self.uprev = y0(self.x)
        self.ucurr = y0(self.x)
        self.alpha = T / m * self.dt ** 2 / self.dx **2
        
        
    def step(self):
        """Steps the system forward in time."""
        unext = np.zeros_like(self.ucurr)
        unext[1: -1] = 2 * self.ucurr[1:-1] - self.uprev[1:-1] + \
                        self.alpha * (self.ucurr[2:] - 2 * self.ucurr[1:-1] + self.ucurr[:-2])
        self.ucurr, self.uprev = unext, self.ucurr
        
    def steps(self, n):
        for _ in range(n):
            self.step()

To test our implementation, we will consider a triangular initial string with the properties of a low E string of a guitar (such as those that you can find on the back of electric guitar strings, as here).

In [2]:
L = 0.6477 # m
T = 7.93 * 9.81 # kg * m.s^-2
f0 = 82.4 # Hz
c = 2 * L * f0
m = T / c**2

triangle = lambda x: (x * 2 / L) * (x < L/2) + (-x * 2 / L + 2) * (x >= L/2) 

td_solver = TimeDomainSolver(L, 200, T, m, triangle)

Let's visualize the initial position of the string:

In [3]:
import holoviews as hv
hv.extension('bokeh')
hv.opts.defaults(hv.opts.Curve(width=500))