8.1

Consider the 1D wave equation

$$ \frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2} $$

8.1 - A-E

8.1 - F

Numerically solve the wave equation for the evolution from an initial condition with u = 0 except for one nonzero node, and verify the stability criterion.

8.1 - G

If the equation is replaced by

$$ \frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2} + \gamma \frac{\partial}{\partial t} \frac{\partial^2 u}{\partial x^2} $$

assume that

$$ u(x, t) = Ae^{i(kx−\omega t)} $$

and find a relationship between $k$ and $\omega$, and simplify it for small $\gamma$. Comment on the relationship to the preceeding question.

8.1 - H

Repeat the numerical solution of the wave equation with the same initial conditions, but include the damping term.

const M = 1000 # time steps
const N = 100 # space steps
state = zeros(M, N)
# set the state for the first 2 steps
state[[1,2], N/4] = 1

const v = 3.0
const dt = 0.01
const dx = 0.1

# n is the current state
for n in 2:(M-1)
    for j in 2:(N-1)
        state[n+1, j] = 2state[n, j] - state[n-1, j] + v^2 * dt^2 / dx^2 * (state[n, j+1] - 2state[n, j] + state[n, j-1])
    end
end

for i in 1:10
    plot_string(state[2i, :], 10, i)
end

using PyCall
using PyPlot
@pyimport matplotlib.animation as anim

function plot_string(x, max, i)
    subplot(max,1,i)
    plot(x')
    ylim(-1, 1)
    yticks([-1, 0, 1])
end

function anim_string(x)
    fig=figure()
    ims = [plot(x[i,:]', c="black") for i=1:size(x,1)]
    ani = anim.ArtistAnimation(fig, ims, interval=25, blit=true)    
    ani[:save]("string.mp4", extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])
    display("text/html", string("""<video autoplay controls><source src="data:video/x-m4v;base64,""",
                                base64(open(readbytes,"string.mp4")),"""" type="video/mp4"></video>"""))
    return nothing
end

anim_string (generic function with 1 method)

sqrt(3/2)

1.224744871391589

Below is an animation with $v = 1.0$, $\Delta t = 0.01$, and $\Delta x = 0.1$

using the stability criterion

$$ \begin{array} \\ \sqrt{\frac{3}{2}} &\ge \left|\frac{v \Delta t}{\Delta x}\right| \\ 1.225 &\ge \frac{1.0 * 0.01}{0.1} \\ 1.225 &\ge 0.1 \end{array} $$

So this solution should be stable

anim_string(state)

Below is an animation with $v = 0.75$, $\Delta t = 0.05$, and $\Delta x = 0.03$

using the stability criterion

$$ \begin{array} \\ \sqrt{\frac{3}{2}} &\ge \left|\frac{v \Delta t}{\Delta x}\right| \\ 1.225 &\ge \frac{0.75 * 0.05}{0.03} \\ 1.225 &\ngeq 1.25 \end{array} $$

So this solution should not be stable. Interestingly, it seems that the threshold for stability is actually when $\left|\frac{v \Delta t}{\Delta x}\right| = 1$, indicating a possible error in my math.

const M = 200 # time steps
const N = 100 # space steps
state = zeros(M, N)
# set the state for the first 2 steps
state[[1,2], N/4] = 1

const v = 0.75
const dt = 0.05
const dx = 0.03

# n is the current state, n+1 is the one we're trying to calculate
for n in 2:(M-1)
    for j in 2:(N-1)
        state[n+1, j] = 2state[n, j] - state[n-1, j] + v^2 * dt^2 / dx^2 * (state[n, j+1] - 2state[n, j] + state[n, j-1])
    end
end

#for i in 1:10
#    plot_string(state[i, :], 10, i)
#end

anim_string(state)

Warning: redefining constant v

For fun, here we start out with a half-sin as the initial condition (for the first two frames) instead of a single nonzero sample.

const M = 1000 # time steps
const N = 100 # space steps
state = zeros(M, N)
# set the state for the first 2 steps
state[1, 1:20] = sin((0:19) / 19 * pi)
state[2, 1:20] = sin((0:19) / 19 * pi)

const v = 5.0
const dt = 0.01
const dx = 0.1

# n is the current state
for n in 2:(M-1)
    for j in 2:(N-1)
        state[n+1, j] = 2state[n, j] - state[n-1, j] + v^2 * dt^2 / dx^2 * (state[n, j+1] - 2state[n, j] + state[n, j-1])
    end
end

anim_string(state)

Here we leave the 1st frame as all zero, and only the 2nd frame is a half-sin. This gives a large time derivative

const M = 1000 # time steps
const N = 100 # space steps
state = zeros(M, N)
# set the state for the first 2 steps
#state[1, 1:20] = sin((0:19) / 19 * pi)
state[2, 1:20] = sin((0:19) / 19 * pi)

const v = 5.0
const dt = 0.01
const dx = 0.1

# n is the current state
for n in 2:(M-1)
    for j in 2:(N-1)
        state[n+1, j] = 2state[n, j] - state[n-1, j] + v^2 * dt^2 / dx^2 * (state[n, j+1] - 2state[n, j] + state[n, j-1])
    end
end

anim_string(state)

const M = 1000 # time steps
const N = 100 # space steps
u = zeros(M, N)
# set the state for the first 2 steps
#state[1, 1:20] = sin((0:19) / 19 * pi)
u[[1,2], N/4] = 1

const v = 3.0
const dt = 0.01
const dx = 0.1
const gamma = 0.005

# n is the current state
for n in 2:(M-1)
    for j in 2:(N-1)
        u[n+1, j] = 2u[n, j] - u[n-1, j] + dt^2/dx^2 * ((v^2+gamma/dt) * (u[n, j+1] - 2u[n, j] + u[n, j-1]) -
        gamma/dt * (u[n-1,j+1] - 2u[n-1, j] + u[n-1, j-1]))
    end
end

anim_string(u)

Warning: redefining constant gamma