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)+ [v2(Δt)2/ (Δx)2][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 xk = a + kΔx. Pick a time interval Δt. Define variables yk(t) = f(xk,t). Pick a small value of Δt. Use above equation to evolve y1, . . . , yN-1 in time, but set y0(t) = yN(t) = 0. To implement our initial conditions, take yk(0) = f0(xk) and yk(-Δt) = f0(xk) - (Δt)g0(xk).
(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 f0(x) = max(1 - |x - 2|, 0) and g0(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 f0(x) as before but try different functions g0(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 f0(x) = sin (kΠ(x-a)/b-a) and g0(x) = 0. What happens? How fast?
(7) The conditions y0(t) = 0 and yN(t) -= 0 implement Dirichlet boundary conditions. To implement Neumann boundary conditions, take y0(t) = y1(t) and YN(t) = YN-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!