20  Wave Equation

Lecture 20

Note

I encourage you to play around with the juptyer notebook for this lecture - you can copy the code with the this notebook button on the side of this page.

Warning

These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.

Tip

This chapter is mostly based on

• Burden, Faires, and Burden (2015) Chapter 11,
  
• Driscoll and Braun (2018) Chapter 10,
• Lecture 16: Driscoll and Braun (2018) Galerkin method,
  
• Lecture 17: Chapter 14.4 of Suli and Mayers (2012): Theorem 14.6.
• Lecture 19: Olver (2014) Chapter 5

20.1 Recap

Assignment 6: You are looking at simple finite difference approximations to the linear transport equation u_t + c u_x = 0. You will see convergence under an assumption on the Courant number \sigma := c\frac{\tau}{h}.

Last time: Convergence of FTCS and BTCS finite difference methods for solving the heat equation u_t - \alpha^2 u_{xx} = 0. For FTCS, we saw that we have convergence when 0 \leq \mu \leq \frac{1}{2} where \mu := \tau (\frac{\alpha}{h})^2 is the diffusion number (and we get faster convergence when \mu = \frac16). Then, we saw that BTCS is unconditionally stable.

Today, we look at the wave equation:

\begin{align} &\square u := u_{tt} - c^2 u_{xx} = 0 \tag{wave}\\ &\begin{split} &u(x,0) = u_0(x) \\ &u_t(x,0) = v_0(x) \end{split} \tag{IC} \\ &\begin{split} &u(a,t) = g_a(t) \\ &u(b,t) = g_b(t) \end{split} \tag{BC} \end{align}

where c is the wave speed. (Note: because the equation is second order in time, we require two initial conditions).

Exercise 20.1 (Ignoring the initial and boundary conditions) show that u_0(x \pm c t) solve the \text{(wave)}. That is, the wave equation supports waves travelling in both directions.

20.2 Finite differences

We again consider a uniform mesh:

\begin{align} x_i &= a + \frac{b-a}{m} i =: a + ih &&\text{for }i=0,\dots,m\nonumber\\ t_j &= \frac{T}{n} j =: j\tau &&\text{for }j=0,\dots,n \end{align}

and aim to approximate u(x_i,t_j) \approx u_{ij}.

The simplest approximation we can think of is to replace the second derivatives with our usual centred difference:

\begin{align} \frac {u_{i,j-1} - 2 u_{ij} + u_{i,j+1}} {\tau^2} % - c^2 \frac {u_{i-1,j} - 2 u_{ij} + u_{i+1,j}} {h^2} % = 0 \end{align}

which is equivalent to

\begin{align} u_{i,j+1} % = \sigma^2 u_{i-1,j} + 2( 1 - \sigma^2) u_{ij} + \sigma^2 u_{i+1,j} % - u_{i,j-1} \end{align}

where \sigma := c\frac{\tau}{h}.

We can impose the boundary conditions by setting u_{0,j} := g_a(t_j) and u_{m,j} := g_b(t_j).

Exercise 20.2 Write the above as a matrix equation for the numerical solution in the interior mesh points \{ u_{i,j+1} \}_{i=1}^{m-1}.

There is a slight problem here! In order to compute u_{i,j+1} we require the solution at the previous two time steps. In order to start the iteration we therefore need to have \{u_{i,0}\} and \{u_{i,1}\}. As we have done so far, we can set u_{i,0} := u_0(x_i) but we also need to specify \{u_{i,1}\}. We do this by using the other initial condition: u_t(x,0) = v_0(x).

Let u be the exact solution to the wave equation. Then, we have

\begin{align} u(x_i,t_1 ) &= u(x_i,\tau) % = u(x_i,0) + \tau u_t(x_i,0) + \frac{\tau^2}{2} u_{tt}(x_i,0) + O(\tau^3) \nonumber\\ % &= u_0(x_i) + \tau v_0(x_i) + \frac{c^2 \tau^2}{2} u_0''(x_i) + O(\tau^3) \nonumber\\ % &= u_0(x_i) + \tau v_0(x_i) + \tfrac{\sigma^2}{2} \left[ u_0(x_{i-1}) - 2u_{0}(x_i) + u_0(x_{i+1}) \right] \nonumber\\ &\qquad + O(\tau^3) + O(h^2) \end{align}

In particular, we could define \{u_{i,1}\} by

\begin{align} u_{i,1} := \tfrac{\sigma^2}{2} u_{i-1,0} +(1- \sigma^2) u_{i,0} + \tfrac{\sigma^2}{2} u_{i+1,0} + \tau v_0(x_i) \end{align}

We implement this here:

function fdm_wave( m, n, c, u0, v0; a=0., b=1., f=0., g_a = 0., g_b = 0., T=.2 )
    # space discretisation
    h = (b-a)/m
    x = collect(a:h:b)

    # time discretisation
    τ = T/n 
    t = collect(0:τ:T)

    σ = c*τ/h
    println( "σ = $σ" )
    
    u = zeros( m+1, n+1 )

    # BC
    u[1,:] = typeof(g_a) <: Function ? g_a.(t) : fill(g_a,n+1)
    u[m+1,:] = typeof(g_b) <: Function ? g_b.(t) : fill(g_b,n+1)

    # IC
    u[:,1] = u0.(x)
    u[2:m,2] = (σ^2/2)u0.(x[1:m-1]) + (1-σ^2)u0.(x[2:m]) +(σ^2/2)u0.(x[3:m+1]) + τ*v0.(x[2:m])

    for j ∈ 2:n
        for i ∈ 2:m
            source = typeof(f) <: Function ? f(x[i],t[i]) : f
            u[i,j+1] = (σ^2)u[i-1,j] + 2(1-σ^2)u[i,j] + (σ^2)u[i+1,j] - u[i,j-1] + τ^2 * source
        end
    end

    return x,t,u
end

c = 1.

# initial condition
u0(x) = exp( -400(x-.3)^2 )
v0(x) = 0.

# exact solution
g(z) = mod(z+1,2)-1 < 0 ? -u0(-(mod(z+1,2)-1)) : u0(mod(z+1,2)-1)
u_e(x,t) = (1/2)*(g(x-t)+g(x+t))

x,t,u = fdm_wave( 90, 1000, c, u0, v0; T=10. )

anim = @animate for j ∈ 1:5:length(t)
    plot(x, u[:, j]; 
        title = L"Wave equation $u_{tt} - u_{xx} = 0$",
        label = "t = $(t[j])", legend=:topright,
        ylims = [-1,1], m=2)
    plot!( x->u_e(x,t[j]); 
        label="exact solution", linestyle=:dot)
end
gif(anim, "pics/wave-1.gif", fps=5)

σ = 0.9

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\wave-1.gif

x,t,u = fdm_wave( 101, 250, c, u0, v0; T=2.5 )

anim = @animate for j ∈ 1:5:length(t)
    plot(x, u[:, j]; 
        title = L"Wave equation $u_{tt} - u_{xx} = 0$",
        label = "t = $(t[j])", legend=:topright,
        ylims = [-1,1], m=2)
    plot!( x->u_e(x,t[j]); 
        label="exact solution", linestyle=:dot)
end
gif(anim, "pics/wave-2.gif", fps=10)

σ = 1.01

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\wave-2.gif

c = 1.

# initial condition
u0(x) = x*sin(2π*x) #abs(x-.5)<.75 ? (abs(x-.5)<.25 ? 0 : 1/2 ) : 0
v0(x) = 0.

# exact solution
g(z) = mod(z+1,2)-1 < 0 ? -u0(-(mod(z+1,2)-1)) : u0(mod(z+1,2)-1)
u_e(x,t) = (1/2)*(g(x-t)+g(x+t))

x,t,u = fdm_wave( 90, 1000, c, u0, v0; T=10. )

anim = @animate for j ∈ 1:5:length(t)
    plot(x, u[:, j]; 
        title = L"Wave equation $u_{tt} - u_{xx} = 0$",
        label = "t = $(t[j])", legend=:topright,
        ylims = [-1,1], m=2)
    plot!( x->u_e(x,t[j]); 
        label="exact solution", linestyle=:dot)
end
gif(anim, "pics/wave-3.gif", fps=5)

σ = 0.9

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\wave-3.gif

21 Courant–Friedrichs–Lewy (CFL) condition

Exercise 21.1 Use the change of variables \xi(x,t) := x-ct and \eta(x,t) := x+ct to rewrite the wave equation \square u = 0 as u_{\xi\eta} = 0. Hence, prove the d’Alembert formula:

\begin{align} u(x,t) = \frac{u_0(x-ct) + u_0(x+ct)}{2} % +\frac{1}{2c}\int_{x-ct}^{x+ct} v_0(y)\mathrm{d}y. \end{align}

The (exact) domain of dependence at (x,t) is the set of all y for which changing the initial data u_0(y), v_0(y) can affect the solution at (x,t).

Exercise 21.2 Use d’Alembert formula to write down the domain of dpendence for the wave equation.

The numerical domain of dependence at (x_i,t_j) is the set of all x_k for which u_{ij} is affected by u_{k0}.

Exercise 21.3 Consider the finite difference scheme as above: u_{i,0} := u_0(x_i) and

\begin{align} u_{i,1} % &= \tau v_0(x_i) + \tfrac{\sigma^2}{2} u_{i-1,0} +(1- \sigma^2) u_{i,0} + \tfrac{\sigma^2}{2} u_{i+1,0} \nonumber\\ % u_{i,j+1} % &= \sigma^2 u_{i-1,j} + 2( 1 - \sigma^2) u_{ij} + \sigma^2 u_{i+1,j} % - u_{i,j-1} \end{align}

where \sigma := c\frac{\tau}{h}. What is the numerical domain of dependence at (x_i,t_j) of this method?

The Courant–Friedrichs–Lewy (CFL) condition states that, for a numerical method to converge, the exact domain of dependence must be contained in the numerical domain of dependence as h, \tau \to 0.

Exercise 21.4 Under what condition on \sigma does the numerical scheme satisfy the CFL condition?

Exercise 21.5 Does the CLF condition explain the instability of the method that we noticed above?

Exercise 21.6 Complete the von Neumann stability analysis to determine when the numerical scheme is stable.

Hint: Take u_{i,j-1} = e^{\mathrm{i}kx_i}, u_{i,j} =\lambda e^{\mathrm{i}kx_i}, and u_{i,j+1} =\lambda^2 e^{\mathrm{i}kx_i} and show that \lambda satisfies a quadratic equation giving \lambda = \alpha \pm \sqrt{\alpha^2 - 1}. Show that when \sigma satisfies the CFL condition then |\lambda|=1, and, when \sigma violates the CFL condition then one of the \lambda's has absolute value > 1.

In this case the CFL condition distinguishes between stable and unstable constants \sigma. This is not true in general - see A6.

21.1 Semi-discretisation

In this lecture we have treated time in the same way as space. Another approach is the so called method of lines:

Another general approach is semi-discretisation: we again take x_i := a + \frac{b-a}{m} i and h := \frac{b-a}{m} and now approximate \{u(x_i, t)\}_{i=0}^m with the (time-dependent) vector \bm u(t). (This is known as the semi-discretisation step because we have discretised in space but not time).

Then, if one is interested in solving the heat equation (u_t = \alpha^2 u_{xx}), we can replace this with

\begin{align} \bm u'(t) &= \alpha^2 D\bm u(t) \end{align}

where D is a differentiation matrix approximating the second derivative in space. This is a linear, constant coefficient system of ODEs that we can solve directly using techniques from Chapter 1 (also see A2).

Exercise 21.7 Approximate the heat equation \square u := u_{tt} - c^2 u_{xx} = 0 with a linear system of ODEs using the method of lines.

Hint: you could define (i) u_t = y and y_t = c^2 u_{xx} or (ii) u_t = z_x and z_t = c^2 u_{x}. (Here, (ii) gives Maxwell’s equations).

21.2 sine-Gordon

We now look at numerical experiments for the sine-Gordon equation:

\begin{align} u_{tt} - c^2 u_{xx} + \sin u = 0. \end{align}

This is a nonlinear equation.

function fdm_sineGordon(m, n, c, u0, v0; a=0., b=1., g_a=0., g_b=0., T=.2)
    # space discretisation
    h = (b-a)/m
    x = collect(a:h:b)

    # time discretisation
    τ = T/n 
    t = collect(0:τ:T)

    σ = c*τ/h
    println("σ = $σ")

    u = zeros(m+1, n+1)

    # BCs 
    u[1,:] .= typeof(g_a) <: Function ? g_a.(t) : g_a
    u[m+1,:] .= typeof(g_b) <: Function ? g_b.(t) : g_b

    # IC
    u[:,1] = u0.(x)

    # u[•,2]
    for i in 2:m
        # Approximation of u_xx at t=0
        u_xx = (u[i+1,1] - 2u[i,1] + u[i-1,1]) / h^2
        u[i,2] = u[i,1] + τ*v0(x[i]) + (τ^2/2) * (c^2 * u_xx - sin(u[i,1]))
    end

    for j ∈ 2:n
        for i ∈ 2:m
            u[i,j+1] = (σ^2) * u[i-1,j] + 2(1 - σ^2) * u[i,j] + (σ^2) * u[i+1,j] - u[i,j-1] - (τ^2) * sin(u[i,j])
        end
    end

    return x, t, u
end

c = 1.
u0(x) = exp(-500(x-.5)^2)
#4 * atan(exp((x - 0.5) / sqrt(1 - 0.5^2))) 
v0(x) = 0.

x, t, u = fdm_sineGordon(100, 2000, c, u0, v0; T=10.0)

anim = @animate for j ∈ 1:10:length(t)
    plot(x, u[:, j]; 
        title = L"Sine-Gordon Equation: $u_{tt} - c^2 u_{xx} + \sin(u) = 0$",
        label = "t = $(round(t[j], digits=2))", 
        legend = :topright,
        ylims = [-1, 1], 
        lw = 2)
end
gif(anim, "pics/sine-Gordon.gif", fps=15)

σ = 0.5

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\sine-Gordon.gif

c = 1.
u_stat(x) = 4 * atan(exp((x-1/2)/c))
u0(x) = u_stat(x)  #+ sin(x)/10
v0(x) = 0.

x, t, u = fdm_sineGordon(100, 100, c, u0, v0; a=-10, b=10 ,g_a=0., g_b=2π, T=10.0)

anim = @animate for j ∈ 1:10:length(t)
    plot(x, u[:, j]; 
        title = "Sine-Gordon (stationary solution)",
        label = "t = $(round(t[j], digits=2))", 
        legend = :topright,
        #ylims = [-1, 1], 
        lw = 2)
    plot!(u_stat, linestyle=:dot, label="u_0")
end
gif(anim, "pics/sine-Gordon-2.gif", fps=15)

σ = 0.5

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\sine-Gordon-2.gif

c = 1.0
ω = 0.6  
η = sqrt(1 - ω^2)

u0(x) = 0.0
v0(x) = (4 * ω * η * cosh(η * x)) / (η^2 * cosh(η * x)^2 + ω^2)

x, t, u = fdm_sineGordon(300, 4000, c, u0, v0; 
                         a=-20.0, b=20.0, T=40.0)

anim = @animate for j ∈ 1:20:length(t)
    plot(x, u[:, j], 
        title = L"Sine-Gordon breather ($\omega = %$ω$)",
        xlabel = "x", ylabel = "u(x,t)",
        label = "t = $(round(t[j], digits=1))",
        ylims = [-4, 4], lw = 2, color = :blue)
end

gif(anim, "pics/sine-Gordon-3.gif", fps = 20)

σ = 0.075

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\sine-Gordon-3.gif

TipTO-DO

• A6 due a week today
• Prepare your presentations
  
• (Sign-up to office hours if you like)

Previous reading:

• Chapter 1 reading:
  • Burden, Faires, and Burden (2015) sections 5.1 - 5.6,
    
  • Driscoll and Braun (2018) zero stability
• Chapter 2 reading:
  • Recap matrices: sections 7.1 - 7.2 of Burden, Faires, and Burden (2015)
  • Jacobi & Gauss-Seidel: section 7.3 Burden, Faires, and Burden (2015)
    
  • CG: section 7.6 Burden, Faires, and Burden (2015)
• Chapter 3 reading:
  • Sections 11.1 - 11.2 of Burden, Faires, and Burden (2015): shooting methods
  • Lecture 16: Driscoll and Braun (2018) Galerkin method
  • Lecture 17: Chapter 14.4 of Suli and Mayers (2012): Theorem 14.6.
  • Lecture 19: Olver (2014) Chapter 5
  • Driscoll and Braun (2018) 11.1, 11.2, 12.4

Burden, Richard L., Douglas J. Faires, and Annette M. Burden. 2015. Numerical Analysis. 10th ed. CENGAGE Learning.

Driscoll, Tobin A, and Richard J Braun. 2018. Fundamentals of Numerical Computation. https://fncbook.com/.

Olver, Peter J. 2014. Introduction to Partial Differential Equations. Undergraduate Texts in Mathematics. Cham: Springer. https://doi.org/10.1007/978-3-319-02099-0.

Suli, Endre, and David F Mayers. 2012. An Introduction to Numerical Analysis. Cambridge, England: Cambridge University Press.