# Benchmark Kuramoto-Sivashinsky DiffEq

**URL:** <https://discourse.julialang.org/t/benchmark-kuramoto-sivashinsky-diffeq/39196>\
**Category:** General Usage\
**Tags:** diffeq\
**Created:** [May 9, 2020, 8:28pm UTC](https://discourse.julialang.org/t/benchmark-kuramoto-sivashinsky-diffeq/39196 "2020-05-09T20:28:32Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 9, 2020, 8:28pm UTC](https://discourse.julialang.org/t/benchmark-kuramoto-sivashinsky-diffeq/39196/1 "2020-05-09T20:28:32Z")

</div>

Hello,

I have been looking at the notebook [KS FDM Work-Precision Diagrams](https://benchmarks.sciml.ai/html/MOLPDE/ks_fdm_wpd.html) 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`:

```julia
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](https://benchmarks.sciml.ai/html/MOLPDE/ks_fdm_wpd.html)

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 10, 2020, 1:38am UTC](https://discourse.julialang.org/t/benchmark-kuramoto-sivashinsky-diffeq/39196/2 "2020-05-10T01:38:13Z")

</div>

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.