5.1 μs

xxxxxxxxxx
 
begin
    using Plots
    using SparseArrays
end

40.3 s

initCond (generic function with 1 method)

xxxxxxxxxx
 
function initCond(xgrid, sType)
    
    if(sType == 1)
        y = 1/(sqrt(2*pi)) * exp.(xgrid .^ 2 ./ (-2)) #gaussian pluck in center
    
    elseif(sType == 2) 
        y = 0.25*sin.(3*xgrid * (pi/xgrid[end-1])) #standing sinusoidal wave
    end 
    
    return (copy(y), copy(y))
end

213 μs

Wave Equation

$\frac{\partial^{2} y (x, t)}{\partial t^{2}} = c^{2} \frac{\partial^{2} y (x, t)}{\partial x^{2}}$

7.4 μs

Finite Difference

$\frac{y_{i}^{j + 1} - 2 y_{i}^{j} + y_{i}^{j - 1}}{(Δ t)^{2}} = c^{2} \frac{y_{i + 1}^{j} - 2 y_{i}^{j} + y_{i - 1}^{j}}{(Δ x)^{2}}$

Let $α = {(\frac{c Δ t}{Δ x})}^{2}$ ;

$y_{i}^{j + 1} = [α y_{i + 1}^{j} + 2 (1 - α) y_{i}^{j} + α y_{i - 1}^{j}] - y_{i}^{j - 1}$

10.0 μs

Matrix form

${\vec{y}}^{j + 1} = M \cdot {\vec{y}}^{j} - {\vec{y}}^{j - 1}$

$M = (\begin{matrix} 2 (1 - α) & α & 0 & . & . & . & . \\ α & 2 (1 - α) & α & 0 & . & . & . \\ 0 & α & 2 (1 - α) & α & 0 & . & . \\ . & . & . & . & . & . & . \\ . & . & . & . & . & . & . \\ . & . & 0 & α & 2 (1 - α) & α & 0 \\ . & . & . & 0 & α & 2 (1 - α) & α \\ . & . & . & . & 0 & α & 2 (1 - α) \end{matrix})$

18.0 μs

Single string

20.3 μs

n :: Number of points in discrete spatial grid
xi :: Initial co-ordinate value in grid
xf :: Final co-ordinate value in grid
m :: Number of timesteps
dt :: Size of each time step
c :: Speed of wave in string

11.3 μs

waveEquation (generic function with 1 method)

xxxxxxxxxx
 
function waveEquation(n::Int64, xi::Float64, xf::Float64, m::Int64, dt::Float64, c::Float64)
    y = zeros(Float64, (n, m)) #height of string at each point in space (n) and time (m)
    
    dx = (xf-xi)/n #spacing of discrete grid
    xgrid = collect(xi:dx:(xf-dx)) #discretized spatial grid
    
    #CREATING EVOLUTION MATRIX
    α = (c*dt/dx)^2 
    M = SparseArrays.spdiagm(-1 => α*ones(n-1), 0 => 2(1-α)*ones(n), 1 => α*ones(n-1))
    
    #INITIAL CONDITIONS
    y[:, 1], y[:, 2] = initCond(xgrid, 1)
    
    #TIME EVOLUTION LOOP
    for j in 2:(m-1)
        y[:, j+1] = M * y[:, j] - y[:, j-1] 
        
        #SOURCE TERMS
        #y[1, j+1] = 0.1*sin(5*j*dt) #sinusoidal driving source
        #y[1, j+1] = 0.125 * exp((j*dt - 2)^2/ (-2*0.1)) #single gaussian pulse
        
        #BOUNDARY CONDITION
        #y[n, j+1] = y[n-1, j+1] #open end
    end
    
    #GIF CREATION
    anim = @gif for j in 1:5:m
       plot(xgrid, y[:, j], ylim = (-0.3, 0.3)) 
    end
    
    return anim
end

197 μs

This code was written hastily, so it does not elegantly accomodate initial/boundary conditions and source terms.

You can however uncomment the specific lines above and run the function call below manually to switch between them.

29.7 μs

xxxxxxxxxx
 
@time waveEquation(200, -10.0, 10.0, 1000, 0.01, 3.0)

33.5 s

Composite string

60.4 μs

n :: Number of points in discrete spatial grid
xi :: Initial co-ordinate value in grid
xf :: Final co-ordinate value in grid
m :: Number of timesteps
dt :: Size of each time step
c1 :: Speed of wave in string 1
c2 :: Speed of wave in string 2

6.4 μs

waveEquationComposite (generic function with 1 method)

xxxxxxxxxx
 
function waveEquationComposite(n::Int64, xi::Float64, xf::Float64, m::Int64, dt::Float64, c1::Float64, c2::Float64)
    y = zeros(Float64, (n, m)) #height of string at each point in space (n) and time (m)
    
    dx = (xf-xi)/n #spacing of discrete grid
    xgrid = collect(xi:dx:(xf-dx))  #discretized spatial grid
    
    #CREATING EVOLUTION MATRIX
    α = (c1*dt/dx)^2
    β = (c2*dt/dx)^2
    
    nmid = n÷2
    M1 = SparseArrays.spdiagm(-1 => α*ones(nmid-1), 0 => 2(1-α)*ones(nmid), 1 => α*ones(nmid-1))
    M2 = SparseArrays.spdiagm(-1 => β*ones(nmid-1), 0 => 2(1-β)*ones(nmid), 1 => β*ones(nmid-1))
    M = SparseArrays.blockdiag(M1, M2)
    M[nmid, nmid + 1] = α
    M[nmid + 1, nmid] = β
    
    #INITIAL CONDITIONS
    #y[:, 1], y[:, 2] = initCond(xgrid, 1)
    
    #TIME EVOLUTION LOOP
    for j in 2:(m-1)
        y[:, j+1] = M * y[:, j] - y[:, j-1] 
        
        #SOURCE TERMS
        #y[1, j+1] = 0.1*sin(5*j*dt) #sinusoidal driving source
        y[1, j+1] = 0.125 * exp((j*dt - 2)^2/ (-2*0.1)) #single gaussian pulse
        
        #BOUNDARY CONDITION
        #y[n, j+1] = y[n-1, j+1] #open end
    end
    
    #GIF CREATION
    anim = @gif for j in 1:5:m
        plot(xgrid[1:nmid], y[1:nmid, j], ylim = (-0.3, 0.3))
        plot!(xgrid[nmid:end], y[nmid:end, j], ylim = (-0.3, 0.3), lw = 2)
    end
    
    return anim
end

247 μs

xxxxxxxxxx
 
waveEquationComposite(200, -10.0, 10.0, 1000, 0.01, 3.0, 5.0)

9.9 s