test.jl

function ksintegrateNaive(u, Lx, dt, Nt, nplot)
    Nx = length(u)                  # number of gridpoints
    kx = vcat(0:Nx/2-1, 0, -Nx/2+1:-1)  # integer wavenumbers: exp(2*pi*i*kx*x/L)
    alpha = 2 * pi * kx / Lx              # real wavenumbers:    exp(i*alpha*x)
    D = 1im * alpha                   # D = d/dx operator in Fourier space
    L = alpha .^ 2 - alpha .^ 4         # linear operator -D^2 - D^4 in Fourier space
    G = -0.5 * D                      # -1/2 D operator in Fourier space
    Nplot = round(Int64, Nt / nplot) + 1  # total number of saved time steps

    x = collect(0:Nx-1) * Lx / Nx
    t = collect(0:Nplot) * dt * nplot
    U = zeros(Nplot, Nx)

    # some convenience variables
    dt2 = dt / 2
    dt32 = 3 * dt / 2
    A = ones(Nx) + dt2 * L
    B = (ones(Nx) - dt2 * L) .^ (-1)

    Nn = G .* fft(u .* u) # -1/2 d/dx(u^2) = -u u_x, collocation calculation
    Nn1 = Nn

    U[1, :] = u # save initial value u to matrix U
    np = 2     # counter for saved data

    # transform data to spectral coeffs 
    u = fft(u)

    # timestepping loop
    for n = 0:Nt-1
        Nn1 = Nn                       # shift N^{n-1} <- N^n
        Nn = G .* fft(real(ifft(u)) .^ 2) # compute N^n = -u u_x

        u = B .* (A .* u + dt32 * Nn - dt2 * Nn1) # CNAB2 formula

        if mod(n, nplot) == 0
            U[np, :] = real(ifft(u))
            np += 1
        end
    end
    U, x, t
end

using FFTW, Gen
Lx = 128
Nx = 1024
dt = 1 / 16
nplot = 8
Nt = 1600

x = Lx * (0:Nx-1) / Nx
u = cos.(x) .+ 0.1 .* cos.(x / 16) .* (1 .+ 2 * sin.(x / 16))

U, x, t = ksintegrateNaive(u, Lx, dt, Nt, nplot)
;

using PyPlot
pcolor(x, t[1:end-1], U)
xlim(x[1], x[end])
ylim(t[1], t[end])
xlabel("x")
ylabel("t")
title("Kuramoto-Sivashinsky dynamics")