# 17 PDEs: waves in 2D 

## 2D wave equation

The wave equation in 2D can be derived in the same fashion as for the 1D case and reads

$$
\left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}\right)u(x, y, t) - \frac{1}{c^2} \frac{\partial^2 u}{\partial t^2} = 0
$$

or in general

$$
\boldsymbol{\nabla}^2 u(\mathbf{x}, t) - \frac{1}{c^2} \frac{\partial^2}{\partial t^2}u(\mathbf{x}, t) = 0.
$$

For the case of a string or a membrane with areal density $\rho$ and constant tension $T$, the wave velocity is

$$
c = \sqrt{\frac{T}{\rho}}.
$$

## Numerical solution: Finite difference discretization 

As before, we replace the derivatives with their finite difference approximations, hence the name *finite difference method* (FDM) for this approach.

\begin{align}
\left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}\right)u(x, y, t) &\approx
 \frac{u(x+\Delta x, y, t) - 2u(x, y, t) + u(x-\Delta x, y, t)}{\Delta x^2} \notag\\
 &+ \frac{u(x, y + \Delta y, t) - 2u(x, y, t) + u(x, y-\Delta y, t)}{\Delta y^2}
\end{align}

Setting $\Delta x = \Delta y = \delta$ (same spacing in $x$ and $y$) and $u_{ijn} := u(i\delta, j\delta, n\Delta t)$ we obtain the discretized 2D wave equation

$$
u_{i,j,n+1} = 2(1 - 2\beta^2)u_{i,j,n} - u_{i,j,n-1} + \beta^2(u_{i+1,j,n} + u_{i-1,j,n} + u_{i,j+1,n} + u_{i,j-1,n})\\
\beta := \frac{c}{\delta / \Delta t}
$$



## Implementation
Packaged as a function.


```python
import matplotlib.pyplot as plt
import numpy as np

plt.style.use('ggplot')
```


```python
def wave2D(L=0.5, N=50, Dt=1e-4, Nt=100, step=1,
           rho=1.5e-3, tension=5.):

    Lx = Ly = L
    Nx = Ny = N
    delta = L/N

    c = np.sqrt(tension/rho)

    beta = c*Dt/delta
    beta2 = beta**2

    print("c = {0} m/s".format(c))
    print("delta = {0} m,  Dt = {1} s, delta/Dt = {2} m/s".format(delta, Dt, delta/Dt))
    print("beta = {}".format(beta))
    print("beta * sqrt(2) = {}".format(beta * np.sqrt(2)))


    x = np.linspace(0, Lx, Nx+1)  # need N+1!
    y = np.linspace(0, Ly, Ny+1)
    #note: we use u[i*delta, j*delta] = u(x, y)
    X, Y = np.meshgrid(x, y)
   
    def gaussian(x, u0=0.05, x0=None, sigma=0.1*L):
        x0 = np.mean(x) if x0 is None else x0
        return u0/np.sqrt(2*np.pi*sigma**2) * np.exp(-(x-x0)**2 / (2*sigma**2))

    # u[i*delta, j*delta] = u(x, y)
    # displacements at j-1, j, j+1
    u0 = np.zeros((Nx+1, Ny+1))
    u1 = np.zeros_like(u0)
    u2 = np.zeros_like(u0)

    # save array
    u_t = np.zeros((int(np.ceil(Nt/step)), Nx+1, Ny+1))

    # use a boolean mask to define regions
    inside = np.zeros_like(u0, dtype=np.bool)
    inside[1:-1, 1:-1] = True
    boundary = np.logical_not(inside)
    
    # boundary conditions
    def set_boundary_conditions(u):
        u[boundary] = 0
  
    set_boundary_conditions(u0)
    set_boundary_conditions(u1)
    set_boundary_conditions(u2)

    # initial conditions: velocity 0, i.e. no difference between y0 and y1

    u0[inside] = u1[inside] = (gaussian(X)*gaussian(Y))[inside]

    # save initial
    t_index = 0
    u_t[t_index, :, :] = u0
    if step == 1:
        t_index += 1
        u_t[t_index, :, :] = u1

    for nt in range(2, Nt):
        u2[1:-1, 1:-1] = 2*(1-2*beta2)*u1[1:-1, 1:-1] - u0[1:-1, 1:-1] + beta2*(u1[2:, 1:-1] 
                                            + u1[:-2, 1:-1] + u1[1:-1, 2:] + u1[1:-1, :-2])
        u0[:, :], u1[:, :] = u1, u2

        if nt % step == 0 or nt == Nt-1:
            t_index += 1
            u_t[t_index, :, :] = u2       
            print("Iteration {0:5d}".format(nt), end="\r")
    else:
        print("Completed {0:5d} iterations: t={1} s".format(nt, nt*Dt))
        
    return u_t, X, Y, delta, Dt, step

```

In 2D, the Courant criterion for stability is

$$
\beta = \frac{c}{\delta/\Delta t} ≤ \frac{1}{\sqrt{2}}
$$

or 

$$
\sqrt{2} \beta ≤ 1.
$$

If the solution becomes unstable, check $\sqrt{2}\beta$!


```python
u_t, X, Y, delta, Dt, step = wave2D(Nt=1000)
```

    c = 57.735026918962575 m/s
    delta = 0.01 m,  Dt = 0.0001 s, delta/Dt = 100.0 m/s
    beta = 0.5773502691896257
    beta * sqrt(2) = 0.816496580927726
    Completed   999 iterations: t=0.0999 s



```python
times = Dt * np.arange(0, len(u_t))
```

Look at the minimum and maximum values for each frame:


```python
fig = plt.figure(figsize=(6,3))
ax = fig.add_subplot(111)
ax.plot(times, u_t.min(axis=(1, 2)),
         times, u_t.max(axis=(1, 2)))
ax.set_xlabel("time (s)")
fig.tight_layout()
```


    <IPython.core.display.Javascript object>






The min/max values do not appear to grow, indicating stable integration.

#### 2D Animation
For 2D animation to work in a Jupyter notebook, use


```python
%matplotlib notebook
```

If no animations are visible, restart kernel and execute the `%matplotlib notebook` cell as the very first one in the notebook.

We use `matplotlib.animation` to look at movies of our solution:


```python
import matplotlib.animation as animation
```

The `update_wave()` function simply re-draws our image for every `frame`.


```python
u_limits = 0.5*u_t.min(), 0.5*u_t.max()

fig1 = plt.figure(figsize=(5,5))
ax = fig1.add_subplot(111)
ax.set_aspect(1)

def update_wave(frame, data):
    global ax, Dt, u_limits
    ax.clear()
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.contour(X, Y, data[frame], 30, cmap=plt.cm.viridis, norm=plt.Normalize(*u_limits))
    ax.contourf(X, Y, data[frame], 30, cmap=plt.cm.viridis, norm=plt.Normalize(*u_limits), alpha=0.8)
    ax.text(0.1, 0.9, "t = {0:g} s".format(frame*Dt), transform=ax.transAxes)

wave_anim = animation.FuncAnimation(fig1, update_wave, len(u_t), fargs=(u_t, ), 
                                    interval=5, blit=False, repeat_delay=1000)

```


    <IPython.core.display.Javascript object>






If you have `ffmpeg` installed then you can export to a MP4 movie file:


```python
wave_anim.save("square_drum_contours.mp4", fps=30, dpi=300)
```

#### 3D animation




```python
from mpl_toolkits.mplot3d import Axes3D
```


```python
def plot_wireframe(X, Y, u, limits=None, ax=None):
    if ax is None:
        fig = plt.figure()
        ax = fig.add_subplot(111, projection="3d")
    ax.plot_wireframe(X, Y, u)
    ax.set_xlabel(r"position $x$ (m)")
    ax.set_ylabel(r"position $y$ (m)")
    ax.set_zlabel(r"displacement $u$ (m)")
    ax.set_zlim(limits)
    ax.figure.tight_layout()
    return ax
```


```python
plot_wireframe(X, Y, u_t[0])
```


    <IPython.core.display.Javascript object>









    <matplotlib.axes._subplots.Axes3DSubplot at 0x1147fb588>




```python
u_limits = u_t.min(), u_t.max()

fig2 = plt.figure(figsize=(7,7))
ax = fig2.add_subplot(111, projection="3d")
ax.set_aspect(1)

def update_wire(frame, data):
    global ax, u_limits, X, Y
    ax.clear()
    plot_wireframe(X, Y, data[frame], limits=u_limits, ax=ax)
    ax.text(0.1, 0.9, 0.9, "t = {0:g} s".format(frame*Dt), transform=ax.transAxes)

wire_anim = animation.FuncAnimation(fig2, update_wire, len(u_t), fargs=(u_t,), 
                                    interval=5, blit=True, repeat_delay=1000)

```


    <IPython.core.display.Javascript object>







```python
wire_anim.save("square_drum_wire.mp4", fps=30, dpi=300)
```

### Demystifying `meshgrid`
Testing with the initial condition, an isotropic 2D [Gaussian](http://mathworld.wolfram.com/GaussianFunction.html)

$$
g_2(x, y) = \frac{1}{2\pi\sigma^2} \exp\left[-\frac{(x-x_0)^2 + (y-y_0)^2}{2\sigma^2}\right]
$$

Note that the 2D Gaussian can be written as a product of two 1D Gaussians

\begin{align}
g_2(x, y) &= \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left[-\frac{(x-x_0)^2}{2\sigma^2}\right] \cdot 
          \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left[-\frac{(y-y_0)^2}{2\sigma^2}\right]\\
        &= g_1(x) \cdot g_1(y)
\end{align}

Plot $g_2(x, y)$ by (1) constructing it from the product of $g_1$s and (2) from the explicit formula. Use [`np.meshgrid`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.meshgrid.html) to evaluate the function on a grid.

Use an *asymmetric* grid ($L_x=0.5, N_x=50$ and $L_y=1, N_y=100$) to clearly see the shape of the underlying arrays.


```python
L = 0.5
Lx, Ly = 0.5, 1
Nx, Ny = 50, 100

x = np.linspace(0, Lx, Nx+1)  # need N+1!
y = np.linspace(0, Ly, Ny+1)
X, Y = np.meshgrid(x, y)
    
def gaussian2D(x, y, u0=0.05, x0=None, y0=None, sigma=0.1*L):
    x0 = np.mean(x) if x0 is None else x0
    y0 = np.mean(y) if y0 is None else y0
    return u0/(2*np.pi*sigma) * np.exp(-((x-x0)**2 + (y-y0)**2)/(2*sigma**2))

def gaussian(x, u0=0.05, x0=None, sigma=0.1*L):
    x0 = np.mean(x) if x0 is None else x0
    return u0/np.sqrt(2*np.pi*sigma) * np.exp(-(x-x0)**2 / (2*sigma**2))
```

#### Meshgrid 


```python
print(x.shape, y.shape)
print(X.shape, Y.shape)
print(X)
print(Y)
```

    (51,) (101,)
    (101, 51) (101, 51)
    [[ 0.    0.01  0.02 ...,  0.48  0.49  0.5 ]
     [ 0.    0.01  0.02 ...,  0.48  0.49  0.5 ]
     [ 0.    0.01  0.02 ...,  0.48  0.49  0.5 ]
     ..., 
     [ 0.    0.01  0.02 ...,  0.48  0.49  0.5 ]
     [ 0.    0.01  0.02 ...,  0.48  0.49  0.5 ]
     [ 0.    0.01  0.02 ...,  0.48  0.49  0.5 ]]
    [[ 0.    0.    0.   ...,  0.    0.    0.  ]
     [ 0.01  0.01  0.01 ...,  0.01  0.01  0.01]
     [ 0.02  0.02  0.02 ...,  0.02  0.02  0.02]
     ..., 
     [ 0.98  0.98  0.98 ...,  0.98  0.98  0.98]
     [ 0.99  0.99  0.99 ...,  0.99  0.99  0.99]
     [ 1.    1.    1.   ...,  1.    1.    1.  ]]


#### Plot 1D Gaussians
Plot the 1D Gaussians using the meshgrid components.


```python
gaussian(X).shape
```




    (101, 51)




```python
%matplotlib notebook
```


```python
fig = plt.figure()
ax = fig.add_subplot(111)
ax.contourf(X, Y, gaussian(X, sigma=0.1), 30, cmap=plt.cm.viridis)
plt.show()
```


    <IPython.core.display.Javascript object>







```python
gaussian(Y).shape
```




    (101, 51)




```python
ax = plt.subplot(131)
ax.contourf(X, Y, gaussian(X, sigma=0.1), 30, cmap=plt.cm.viridis)
ax.set_aspect(1)
ax = plt.subplot(132)
ax.contourf(X, Y, gaussian(Y, sigma=0.1), 30, cmap=plt.cm.viridis)
ax.set_aspect(1)
ax = plt.subplot(133)
ax.contourf(X, Y, gaussian(X, sigma=0.1) * gaussian(Y, sigma=0.1), 30, cmap=plt.cm.viridis)
ax.set_aspect(1)
```

#### 2D gaussian 


```python
gaussian2D(X, Y).shape
```




    (101, 51)




```python
ax = plt.subplot(111)
ax.contourf(X, Y, gaussian2D(X, Y, sigma=0.1), 30, cmap=plt.cm.viridis)
ax.set_aspect(1)
```

d9b13e75e31883eb7353c5f075abc30823a94711
ipynb
Jupyter Notebook
17_PDEs_waves/17_PDEs_2D_waves.ipynb
nachrisman/PHY494
bac0dd5a7fe6f59f9e2ccaee56ebafcb7d97e2e7
CC-BY-4.0
17_PDEs_waves/17_PDEs_2D_waves.ipynb
nachrisman/PHY494
bac0dd5a7fe6f59f9e2ccaee56ebafcb7d97e2e7
CC-BY-4.0
17_PDEs_waves/17_PDEs_2D_waves.ipynb
nachrisman/PHY494
bac0dd5a7fe6f59f9e2ccaee56ebafcb7d97e2e7
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES