Benchmark Kuramoto-Sivashinsky DiffEq

General Usage
1

Hello,

I have been looking at the notebook KS FDM Work-Precision Diagrams that solves the Kuramot-Sivashinsky (KS) problem with finite differences.

The 1D equation is:

\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \frac{\partial^2 u}{\partial x^2} + \frac{\partial^4 u}{\partial x^4} = 0.

I am not sure about certain coefficients in lin_term:

DiffEqArrayOperator(-0.0004*((1/dx^2) * diagm(-1 => du2, 0 => d2, 1 => du2)
                        +(1/dx^4) * diagm(-2 => duu4, -1 => du4, 0 => d4, 1 => du4, 2 => duu4)))

Why is there a -0.0004 coefficient?
Also, I would expect to have a minus sign in '(1/dx^4) * diagm(-2 => duu4, -1 => du4, 0 => d4, 1 => du4, 2 => duu4))``` since this term got moves to the right hand side .

The full code from KS FDM Work-Precision Diagrams

using ApproxFun, OrdinaryDiffEq, Sundials
using DiffEqDevTools
using LinearAlgebra
using Plots; gr()

# Define the linear and nonlinear terms
function lin_term(N)
    #which is -(D2+D4)
    dx = 1/(N + 1)
    d2 = (-2) * ones(N) # main diagonal
    du2 = ones(N - 1) # off diagonal

    d4 = 6 * ones(N) # main diagonal
    du4 = (-4) * ones(N - 1) # off diagonal
    duu4 = ones(N - 2)
    DiffEqArrayOperator(-0.0004*((1/dx^2) * diagm(-1 => du2, 0 => d2, 1 => du2)
                        +(1/dx^4) * diagm(-2 => duu4, -1 => du4, 0 => d4, 1 => du4, 2 => duu4)))
end

function nl_term(N)
    dx = 1/(N + 1)
    du = ones(N - 1) # super diagonal
    dl = -ones(N - 1) # lower diagonal
    D = (-0.2/(4*dx)) * diagm(-1 => dl, 1 => du)

    tmp = zeros(N)
    function (du,u,p,t)
        @. tmp = u^2
        mul!(du, D, tmp)
    end
end

# Construct the problem
function ks(N)
    f1 = lin_term(N)
    f2 = nl_term(N)
    dx = 1 / (N + 1)
    xs = (1:N) * dx

    μ0 = 0.3; σ0 = 0.05
    f0 = x -> 0.6*exp(-(x - μ0)^2 / (2 * σ0^2))
    u0 = f0.(xs)
    prob = SplitODEProblem(f1, f2, u0, (0.0, 1.0))
    xs, prob
end;


xs, prob = ks(100)
sol = solve(prob, RadauIIA5(autodiff=false); abstol=1e-14, reltol=1e-14)
test_sol = TestSolution(sol);

tslices = [0.0 0.25 0.50 0.75 1.]
ys = hcat((sol(t) for t in tslices)...)
labels = ["t = $t" for t in tslices]
plot(xs, ys, label=labels)
2

The equation can have arbitrary coefficients of course, but we should’ve documented this. Open an issue and we’ll get it fixed up. Technically you can say the 4th order derivative term has a -1 coefficient, but we probably should have it be -. In general, a lot of these PDE benchmarks need to be updated, and that just takes a bit of time. The student who was working on them was just there for a summer, but left a pretty good template for someone to really finish.

2 Likes