The problems that we've dealt with thus far all model natural diffusion processes: heat transport, molecular diffusion, etc. Another interesting physical phenomenon is that of wave propagation. In 1 spatial dimension the wave equation is \begin{gather*} u_{tt} = \alpha^2 u_{xx} \end{gather*} where \(\alpha\) is the stiffness of the wave. With Dirichlet boundary conditions we can think of this as the behaviour of a guitar string after it has been plucked.

Let's write MATLAB code to numerically solve this problem:

Consider \(u_{tt} = 2 u_{xx}\) in \(x \in (0,1)\) with \(u(0,x) = x(1-x)\), \(u_t(0,x) = 0\), and \(u(t,0) = u(t,1) = 0\) and \(\alpha = 5\). We can discretize the derivatives as \begin{align*} u_{tt}(t_{n+1},x_j) \amp \approx \frac{U_j^{n+1} - 2U_j^n + U_j^{n-1}}{\Delta t^2}\\ u_{xx}(t_n,x_j) \amp \approx \frac{U_{j+1}^n - 2U_j^n + U_{j-1}^n}{\Delta x^2} \end{align*}

There is a natural stability question waiting to be asked about the discretization of the 1D wave equation. Ask and answer this question.

Write code to solve the 2D Dirichlet problem \begin{equation*} \text{ Solve: } u_{tt} = \alpha^2 \nabla \cdot \nabla u \text{ in } x \in \Omega \text{ with } u(0,x,y) = \sin(2\pi (x-0.5))\sin(2\pi(y-0.5)) \end{equation*} subject to homogeneous Dirichlet boundary conditions.

For simplicity we suggest that you take \(\Delta x = \Delta y\). You should also be very careful of the stability conditions for the heat equation.

A traveling wave can be modeled by the PDE \begin{equation*} u_t + a u_x = 0 \end{equation*} where \(a\) is the speed of the wave propogation and \(u(t,x)\) is the height of the wave. Write a numerical solve for the traveling wave problem on \(x \in (0,\infty)\) with initial condition \(u(0,x) = \text{ exp } \left( -\frac{(x-1)^2}{0.1} \right)\), boundary condition \(u(t,0) = 0\), and \(a=1\).

Solve this problem with and Euler-type time step and

1. centered differences in space, and

2. upwind differences in space.

where the “wind” is from left to right. Plot both solutions on top of each other. What do you notice about the behavior of the solutions? Neither of these solutions should actually give a traveling wave, but that is what is expected out of the solution.

Note about the analytic solution: If \(\eta(x)\) is the initial condition for \(u_t + a u_x = 0\) then \(u(t,x) = \eta(x-at)\) is the analytic solution to the PDE. You should check this using the chain rule.

Three ways to fix the issues seen in the previous problem are the “Leapfrog” scheme, the “Lax-Friedrichs” scheme, and the “Lax-Wendroff” scheme: \begin{align*} \text{ Leapfrog: } \amp \frac{U_j^{n+1} - U_j^{n-1}}{2\Delta t} = -a \frac{U_{j+1}^n - U_{j-1}^n}{2\Delta x}\\ \text{ Lax-Friedrichs: } \amp \frac{U_j^{n+1} - \frac{1}{2} \left( U_{j+1}^n + U_{j-1}^n \right)}{\Delta t} = -a \frac{U_{j+1}^n - U_{j-1}^n}{2\Delta x}\\ \text{ Lax-Wendroff: } \amp U_j^{n+1} = U_j^n - \frac{a \Delta t}{2\Delta x} \left( U_{j+1}^n - U_{j-1}^n \right) + \left( \frac{a^2 \Delta t^2}{2\Delta x^2} \right) \left( U_{j-1}^n - 2 U_j^n + U_{j+1}^n. \right) \end{align*}

Implement all of these schemes and discuss stability and consistency.