Wave Equation: Program 1

Solve the 1-D Wave Equation

\[ \frac{\partial^2 U}{\partial t^2} = v^2\frac{\partial^2 U}{\partial x^2} \]

Boundary Conditions:

\( U(0, t) = U(10, t) = 0 \quad (t > 0) \)

Initial Conditions:

\[ U(x, 0) = e^{-k(x - x_0)^2} \quad \text{for } 0 \leq x \leq 10 \]
\[ \left. \frac{\partial U}{\partial t} \right|_{x, t = 0} = 0 \quad \text{for } 0 \leq x \leq 10 \]

Take \( x_0 = 5.0 \)
Take \( k = 10.0 \)

Plot the numerical solution \( U_n(x, t) \) with respect to \( x \) for any four time instants within \( 0 \leq t \leq 15 \).

import numpy as np
import matplotlib.pyplot as plt

# Constants
L = 10.0          # Length of the rod
x0 = 5.0          # Initial position of the spike
k = 10            # Gaussian decay constant
Tf = 15.0         # Final time
v = 1.0           # Wave speed

# Discretization
x, dx = np.linspace(0, L, 1001, retstep=True)
t, dt = np.linspace(0, Tf, 4001, retstep=True)

# Stability parameters
r = (v * dt) / dx
c = r**2

# Boundary condition initialization
U = np.zeros((len(x), len(t)))
U[0, :] = 0.0           # Boundary condition at x=0
U[-1, :] = 0.0          # Boundary condition at x=L
U[:, 0] = np.exp(-k * (x - x0)**2)  # Initial condition

# First time step calculation
for i in range(1, len(x) - 1):
    U[i, 1] = 0.5 * (c * U[i + 1, 0] + 2 * (1 - c) * U[i, 0] + c * U[i - 1, 0])

# Time-stepping loop
for n in range(1, len(t) - 1):
    for i in range(1, len(x) - 1):
        U[i, n + 1] = c * U[i + 1, n] + 2 * (1 - c) * U[i, n] + c * U[i - 1, n] - U[i, n - 1]

# Plotting results at specified time intervals
for t_index in range(0, len(t), 1100):
    plt.plot(x, U[:, t_index], ls='-.', lw=2, label=f't = {round(t_index * dt, 1)} s')

plt.xlabel('x')
plt.ylabel('U')
plt.legend()
plt.tight_layout()
plt.savefig('plot.png')

Output 1

The 1D Wave Equation

The general form of the 1D wave equation is:

\[ \frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2} \]

1. Variables

\( u(x,t) \) is the displacement (or wave height) at position \( x \) and time \( t \)
\( v \) is the wave speed

2. Discretization

We discretize both space and time:

\( u_i^n \) represents the value of \( u(x,t) \) at the \( i \)-th spatial grid point \( x_i \) and \( n \)-th time step \( t_n \).
The grid spacing in space is \( \Delta x = dx \) and in time is \( \Delta t = dt \).

Spatial Derivatives:

Using a central difference approximation for the second spatial derivative:

\[ \frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{(\Delta x)^2} \]

Time Derivatives:

Similarly, the second-order time derivative can be approximated using a central difference:

\[ \frac{\partial^2 u}{\partial t^2} \approx \frac{u_i^{n+1} - 2u_i^n + u_i^{n-1}}{(\Delta t)^2} \]

3. Recurrence Relation from Wave Equation

From the wave equation:

\[ \small \frac{u_i^{n+1} - 2u_i^n + u_i^{n-1}}{(\Delta t)^2} = v^2 \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{(\Delta x)^2} \]

Rearranging this equation to solve for \( u_i^{n+1} \):

\[ \small u_i^{n+1} = 2u_i^n - u_i^{n-1} + \left(\frac{v^2 (\Delta t)^2}{(\Delta x)^2}\right)\left(u_{i+1}^n - 2u_i^n + u_{i-1}^n\right) \]

4. Initial Conditions and Boundary Conditions

Boundary conditions are set at \( x = 0 \) and \( x = L \) (i.e., \( u_0^n = u_{N-1}^n = 0 \), for all \( n \)).

The initial condition is set as a Gaussian spike at \( t = 0 \):

\[ u(x, 0) = \exp(-k(x - x_0)^2) \]

5. First Time Step Approximation

To start with, we have for j = 0

\[ \small u_i^n = 2u_i^0 - u_i^{-1} + \left(\frac{v^2 (\Delta t)^2}{(\Delta x)^2}\right)\left(u_{i+1}^0 - 2u_i^0 + u_{i-1}^0\right) \]

As the initial condition \[ \left. \frac{\partial{u}}{\partial{t}}\right|_{t=0}=0 \], we can write \[ \frac{u_i^1 - u_i^{-1}}{2 \Delta t} = 0,\quad i.e. \quad u_i^1 = u_i^{-1} \]

At \( n = 0 \), the first time step is calculated differently using:

\[ u_i^1 = 0.5 \cdot \left( c \cdot (u_{i+1}^0 + u_{i-1}^0) + 2(1 - c) \cdot u_i^0 \right) \]

6. Final Recurrence Relation

The final recurrence relation implemented in the time-stepping loop for (n > 0)is:

\[ \small u_i^{n+1} = c \cdot u_{i+1}^n + 2(1 - c) \cdot u_i^n + c \cdot u_{i-1}^n - u_i^{n-1} \]

Where:

\[ c = \left(\frac{v \Delta t}{\Delta x}\right)^2 \]

Partial Differential Equations