Exploration: Discretizing the Wave Equation

A computer cannot store an entire function in its memory. Typically it stores the values of the function at a finite number of points. It approximates derivatives by differences of functions. In particular,

Plugging this into equation gives a discretized version of the wave equation:

F(x, t + Δt) = 2f (x, t) - f (x, t - Δt)+ [v²(Δt)²/ (Δx)²][f(x + Δx, t) - 2f(x, t) + f(x- Δx, t)]

We will use above equation to solve the wave equation numerically. Take an interval [a, b] and a large integer N, let Δx = (b - a)/N, and let x_k = a + kΔx. Pick a time interval Δt. Define variables y_k(t) = f(x_k,t). Pick a small value of Δt. Use above equation to evolve y₁, . . . , yN-1 in time, but set y₀(t) = y_N(t) = 0. To implement our initial conditions, take y_k(0) = f₀(x_k) and y_k(-Δt) = f₀(x_k) - (Δt)g₀(x_k).

(1) Pick values of a, b, v, N, and Δt, and write a program to implement this process. Be sure to leave a, b, v, N, and Δt as parameters that can be changed easily.

(2) Run this program with a = 0, b = 10, v = 1, N = 100, and Δt = 0.01 with initial conditions f₀(x) = max(1 - |x - 2|, 0) and g₀(x) = 0. Plot f (x, t) for t ranging from 0 to 20 at intervals of 0.5.

(3) Vary the sizes of Δt and Δx, always keeping vΔt  <<  Δx. How small do Δx and Δt need to be to adequately model the continuous wave equation?

(4) Now take Δx <  vΔt. What happens? Why? If your program keeps crashing, take Δx only slightly less than vΔt.

Viewing the discretized wave equation as a second-order difference equation in N-1 variables, one of the modes crosses the line from neutrally stable to unstable when Δx hits vΔt. As Δx gets smaller and smaller, this mode becomes more and more unstable, and other modes become unstable as well. This is to be expected, since the discretized wave equation can only carry a pulse a distance Δx in time Δt. If Δx <  vΔt, the discretized wave equation cannot support traveling waves.

(5) Going back to the setup of step 2, take f₀(x) as before but try different functions g₀(x). Can you find a function for which the wave is purely a forward traveling wave at time t = 0?

(6) For various integers k, take f₀(x) = sin (kΠ(x-a)/b-a) and g₀(x) = 0. What happens? How fast?

(7) The conditions y₀(t) = 0 and y_N(t) -= 0 implement Dirichlet boundary conditions. To implement Neumann boundary conditions, take y₀(t) = y₁(t) and Y_N(t) = Y_N-1(t). Repeat step 2 with these boundary conditions.

(8) What happens if you take Dirichlet boundary conditions at x = a and Neumann boundary conditions at x = b? Try it!