# Benjamin-Bona-Mahony equation


```julia
include("setup.jl")
```

# Periodic boundary conditions


```julia
function relaxation_functional(γ, unew, uold, param)
    @unpack D, D2, tmp1, tmp2 = param
    
    @. tmp2 = (1-γ)*uold + γ*unew
    mul!(tmp1, D2, tmp2)
    @. tmp1 = tmp2^2 - tmp2 * tmp1
    energy = integrate(tmp1, D)
end
function relaxation_functional(u, param)
    @unpack D, D2, tmp1 = param
    
    mul!(tmp1, D2, u)
    @. tmp1 = u^2 - u * tmp1
    energy = integrate(tmp1, D)
end

function relaxation!(integrator)
    told = integrator.tprev
    uold = integrator.uprev
    tnew = integrator.t
    unew = integrator.u
    
    γ = one(tnew)
    terminate_integration = false
    γlo = one(γ)/2
    γhi = 3*one(γ)/2
    energy_old = relaxation_functional(uold, integrator.p)
    if (relaxation_functional(γlo, unew, uold, integrator.p)-energy_old) * (relaxation_functional(γhi, unew, uold, integrator.p)-energy_old) > 0
        terminate_integration = true
    else
        γ = find_zero(g -> relaxation_functional(g, unew, uold, integrator.p)-energy_old, (γlo, γhi), Roots.AlefeldPotraShi())
    end
    if γ < eps(typeof(γ))
        terminate_integration = true
    end
    
    @. unew = uold + γ * (unew - uold)
    DiffEqBase.set_u!(integrator, unew)
    if !(tnew ≈ top(integrator.opts.tstops))
        tγ = told + γ * (tnew - told)
        DiffEqBase.set_t!(integrator, tγ)
    end
    
    if terminate_integration
        terminate!(integrator)
    end
    
    nothing
end


function save_func_bbm_periodic(u, t, integrator)
    @unpack D, D2, tmp1, usol = integrator.p
    
    mass = integrate(u, D)
    
    mul!(tmp1, D2, u)
    @. tmp1 = u^2 - u * tmp1
    quadratic = integrate(tmp1, D)
    
    @. tmp1 = (u + 1)^3
    cubic = integrate(tmp1, D)
    
    SVector(mass, quadratic, cubic)
end


function bbm_periodic!(du, u, param, t)
    @unpack D, invImD2, tmp1, tmp2 = param
    one_third = one(t) / 3
    
    # conservative semidiscretization
    @. tmp1 = -one_third * u^2
    mul!(tmp2, D, tmp1)
    mul!(tmp1, D, u)
    @. tmp2 += -one_third * u * tmp1 - tmp1
    ldiv!(du, invImD2, tmp2)
    
    nothing
end

function solve_ode_bbm_periodic(usol, D, D2, tspan, alg, tol, dt, adaptive)
    invImD2 = isa(D2, AbstractMatrix) ? factorize(I - D2) : I - D2
    
    x = grid(D)
    u0 = usol.(tspan[1], x, x[1], -x[1])
    tmp1 = similar(u0); tmp2 = similar(tmp1)
    param = (D=D, D2=D2, invImD2=invImD2, tmp1=tmp1, tmp2=tmp2, usol=usol)

    ode = ODEProblem(bbm_periodic!, u0, tspan, param)
    
    saveat = range(tspan..., length=100)
    saved_values_baseline = SavedValues(eltype(D), SVector{3,eltype(D)})
    saving_baseline = SavingCallback(save_func_bbm_periodic, saved_values_baseline, saveat=saveat)
    saved_values_relaxation = SavedValues(eltype(D), SVector{3,eltype(D)})
    saving_relaxation = SavingCallback(save_func_bbm_periodic, saved_values_relaxation, saveat=saveat)
    relaxation = DiscreteCallback((u,t,integrator) -> true, relaxation!, save_positions=(false,true))
    cb_baseline = CallbackSet(saving_baseline)
    cb_relaxation = CallbackSet(relaxation, saving_relaxation)
    
    sol_relaxation = solve(ode, alg, abstol=tol, reltol=tol, dt=dt, adaptive=adaptive, save_everystep=false,
        callback=cb_relaxation, tstops=saveat)
    sol_baseline = solve(ode, alg, abstol=tol, reltol=tol, dt=dt, adaptive=adaptive, save_everystep=false,
        callback=cb_baseline, tstops=saveat)

    unum_baseline   = sol_baseline[end]
    unum_relaxation = sol_relaxation[end]
    uana = usol.(tspan[end], x, x[1], -x[1])
    @printf("Error in u (baseline):   %.3e\n", integrate(u->u^2, unum_baseline - uana, D) |> sqrt)
    @printf("Error in u (relaxation): %.3e\n", integrate(u->u^2, unum_relaxation - uana, D) |> sqrt)
    @printf("Difference of baseline and relaxation in u: %.3e\n", 
        integrate(u->u^2, unum_baseline - unum_relaxation, D) |> sqrt)

    sleep(0.1)
    fig_u, ax = plt.subplots(1, 1)
    plt.plot(x, u0, label=L"u^0")
    plt.plot(x, uana, label=L"$u^\mathrm{ana}$")
    plt.plot(x, unum_baseline, label=L"$u^\mathrm{num}$ (baseline)")
    plt.plot(x, unum_relaxation, label=L"$u^\mathrm{num}$ (relaxation)")
    plt.xlabel(L"x"); plt.ylabel(L"u")
    plt.legend(loc="center left", bbox_to_anchor=(1.0, 0.5));

    t_baseline = saved_values_baseline.t
    t_relaxation = saved_values_relaxation.t
    mass_baseline      = map(x->x[1], saved_values_baseline.saveval)
    mass_relaxation    = map(x->x[1], saved_values_relaxation.saveval)
    quadratic_baseline   = map(x->x[2], saved_values_baseline.saveval)
    quadratic_relaxation = map(x->x[2], saved_values_relaxation.saveval)
    cubic_baseline       = map(x->x[3], saved_values_baseline.saveval)
    cubic_relaxation     = map(x->x[3], saved_values_relaxation.saveval)

    fig_invariants, ax = plt.subplots(1, 1)
    ax.set_yscale("symlog", linthreshy=1.0e-14)
    plt.plot(t_baseline,   mass_baseline   .- mass_baseline[1], label=L"$\int u$ (baseline)",
                color="#E69F00", linestyle="-")
    plt.plot(t_relaxation, mass_relaxation .- mass_relaxation[1], label=L"$\int u$ (relaxation)",
                color="#56B4E9", linestyle="-")
    plt.plot(t_baseline,   quadratic_baseline   .- quadratic_baseline[1], label=L"$\int (u^2 + (\partial_x u)^2)$ (baseline)",
                color="#E69F00", linestyle="--")
    plt.plot(t_relaxation, quadratic_relaxation .- quadratic_relaxation[1], label=L"$\int (u^2 + (\partial_x u)^2)$ (relaxation)",
                color="#56B4E9", linestyle="--")
    plt.plot(t_baseline,   cubic_baseline    .- cubic_baseline[1], label=L"$\int (u + 1)^3$ (baseline)",
                color="#E69F00", linestyle=":")
    plt.plot(t_relaxation, cubic_relaxation  .- cubic_relaxation[1], label=L"$\int (u + 1)^3$ (relaxation)",
                color="#56B4E9", linestyle=":")
    plt.xlabel(L"t"); plt.ylabel("Change of Invariants")
    plt.legend(loc="center left", bbox_to_anchor=(1.0, 0.5))
        
    fig_u, fig_invariants
end

```


```julia
# traveling wave solution
xmin = -90.
xmax = -xmin
get_c() = 1.2
function usol(t, x, xmin, xmax)
    c = get_c()
    A = 3 * (c-1)
    K = 0.5 * sqrt(1 - 1/c)
    x_t = mod(x - c*t - xmin, xmax - xmin) + xmin
    
    A / cosh(K * x_t)^2
end

println("c = ", get_c())
println("xmin = ", xmin)
println("xmax = ", xmax)
@show usol(0., xmin, xmin, xmax)
@show usol(0., xmax, xmin, xmax)
@show N = 2^8
@show dt = 1.0 * (xmax - xmin) / (N * get_c())
flush(stdout)

# tspan = (0., (xmax-xmin)/(3*get_c()))
# tspan = (0.0, (xmax-xmin)/(3*get_c()) + 1*(xmax-xmin)/get_c())
tspan = (0.0, (xmax-xmin)/(3*get_c()) + 10*(xmax-xmin)/get_c())
# tspan = (0.0, (xmax-xmin)/(3*get_c()) + 100*(xmax-xmin)/get_c())
# tspan = (0.0, (xmax-xmin)/(3*get_c()) + 1000*(xmax-xmin)/get_c())

tol = 1.0e-7
adaptive = false
D = fourier_derivative_operator(xmin, xmax, N)
# D = periodic_derivative_operator(1, 6, xmin, xmax, N+1)
# D = periodic_derivative_operator(Holoborodko2008(), 1, 4, xmin, xmax, N+1)
D2 = D^2
# D2 = periodic_derivative_operator(2, 6, xmin, xmax, N+1)
# D2 = periodic_derivative_operator(Holoborodko2008(), 1, 4, xmin, xmax, N+1)^2
# p = 7; mesh = UniformPeriodicMesh1D(xmin, xmax, N÷p)
# D = couple_continuosly(legendre_derivative_operator(-1., 1., p+1), mesh)
# D2 = sparse(couple_continuosly(legendre_second_derivative_operator(-1., 1., p+1), mesh))
# D2 = sparse(D)^2
# p = 7; Dop = legendre_derivative_operator(-1., 1., p+1); mesh = UniformPeriodicMesh1D(xmin, xmax, N÷(p+1))
# D = couple_discontinuosly(Dop, mesh); D₊ = couple_discontinuosly(Dop, mesh, Val(:plus)); D₋ = couple_discontinuosly(Dop, mesh, Val(:minus))
# D2 = sparse(D₊) * sparse(D₋)
# D2 = sparse(D₋) * sparse(D₊)
# D2 = sparse(D)^2

fig_u, fig_invariants = solve_ode_bbm_periodic(usol, D, D2, tspan, RK4(), tol, dt, adaptive)
```

# Convergence study with manufactured solutions


```julia
import SymPy; sp = SymPy

function math_replacements(s)
    s = replace(s, "cos(pi*" => "cospi(")
    s = replace(s, "sin(pi*" => "sinpi(")
end

function usol(t, x)
    exp(t/2) * sinpi(2*(x-t/2))
end

let (t, x) = sp.symbols("t, x", real=true)
    u = usol(t, x)
    
    println("u:")
    sp.diff(u, t) - sp.diff(u, x, 2, t, 1) + sp.diff(u, x) + sp.diff(u^2/2, x) |> 
        sp.simplify |> sp.string |> math_replacements |> println
end
```


```julia
function save_func_bbm_manufactured(u, t, integrator)
    @unpack D, tmp1, usol = integrator.p
    
    x = grid(D)
    tmp1 .= ( usol.(t, x) .- u ).^2
    error_u = integrate(tmp1, D) |> sqrt
    
    error_u
end

function bbm_periodic_manufactured!(du, u, param, t)
    @unpack D, invImD2, tmp1, tmp2 = param
    x = grid(D)
    one_third = one(t) / 3
    
    # conservative semidiscretization
    @. tmp1 = -one_third * u^2
    mul!(tmp2, D, tmp1)
    mul!(tmp1, D, u)
    @. tmp2 += -one_third * u * tmp1 - tmp1
    @. tmp2 += -2*pi^2*(sinpi((t - 2*x)) + 2*pi*cospi((t - 2*x)))*exp(t/2) - exp(t/2)*sinpi((t - 2*x))/2 + pi*exp(t/2)*cospi((t - 2*x)) - pi*exp(t)*sinpi((2*t - 4*x))
    ldiv!(du, invImD2, tmp2)
    
    nothing
end

function errors_bbm_periodic_manufactured(usol, D, D2, tspan, alg, tol, dt, adaptive)
    invImD2 = isa(D2, AbstractMatrix) ? factorize(I - D2) : I - D2
    
    x = grid(D)
    u0 = usol.(tspan[1], x)
    tmp1 = similar(u0); tmp2 = similar(tmp1)
    param = (D=D, D2=D2, invImD2=invImD2, tmp1=tmp1, tmp2=tmp2, usol=usol)

    ode = ODEProblem(bbm_periodic_manufactured!, u0, tspan, param)
    
    saveat = range(tspan..., length=2)
    saved_values_baseline = SavedValues(eltype(D), eltype(D))
    saving_baseline = SavingCallback(save_func_bbm_manufactured, saved_values_baseline, saveat=saveat)
    cb_baseline = CallbackSet(saving_baseline)
    
    sol_baseline = solve(ode, alg, abstol=tol, reltol=tol, dt=dt, adaptive=adaptive, save_everystep=false,
        callback=cb_baseline, tstops=saveat)
    
    error_u = saved_values_baseline.saveval[end]
end

xmin = 0.0
xmax = 1.0
tspan = (0.0, 1.0)

tol = 1.0e-12

# val_N = round.(Int, 2 .^ range(4, 5, length=2))
# val_N = round.(Int, 2 .^ range(4, 8, length=7)) |> evenodd_values
val_N = round.(Int, 2 .^ range(3, 4.5, length=7)) |> evenodd_values
# val_N = 2 .^ (2:7) .+ 1
val_error_u = Float64[]
for N in val_N
#     D = fourier_derivative_operator(xmin, xmax, N)
#     p = 6; D = periodic_derivative_operator(1, p, xmin, xmax, N+1)
#     p = 4; D = periodic_derivative_operator(Holoborodko2008(), 1, p, xmin, xmax, N+1)
#     D2 = D^2
#     D2 = periodic_derivative_operator(2, p, xmin, xmax, N+1)
#     D2 = periodic_derivative_operator(Holoborodko2008(), 1, p, xmin, xmax, N+1)^2
    p = 3; mesh = UniformPeriodicMesh1D(xmin, xmax, N)
    D = couple_continuosly(legendre_derivative_operator(-1., 1., p+1), mesh)
    D2 = sparse(couple_continuosly(legendre_second_derivative_operator(-1., 1., p+1), mesh))
#     D2 = sparse(D)^2
    p = 6; Dop = legendre_derivative_operator(-1., 1., p+1); mesh = UniformPeriodicMesh1D(xmin, xmax, N)
    D = couple_discontinuosly(Dop, mesh); D₊ = couple_discontinuosly(Dop, mesh, Val(:plus)); D₋ = couple_discontinuosly(Dop, mesh, Val(:minus))
    D2 = sparse(D₊) * sparse(D₋)
    D2 = sparse(D₋) * sparse(D₊)
    D2 = sparse(D)^2

    error_u = errors_bbm_periodic_manufactured(usol, D, D2, tspan, Tsit5(), tol, 1/N, true)
    push!(val_error_u, error_u)
end

fig, ax = plt.subplots(1, 1)
ax.set_prop_cycle(marker_cycler)
plt.plot(val_N, val_error_u, label=L"$\| u - u_{\mathrm{ana}} \|$")
c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
plt.plot(val_N, exp(c0) .* val_N.^c1, marker="", linestyle=":", color="gray", label=@sprintf("Order %.2f", -c1))
plt.xscale("log", basex=2)
plt.yscale("symlog", linthreshy=1.0e-13)
plt.legend(loc="center left", bbox_to_anchor=(1.0, 0.5))
ax.locator_params(axis="y", numticks=9)
```

## Convergence study: Plots FD


```julia
function do_stuff(p, val_N)
    val_error = Float64[]
    for N in val_N
        D = periodic_derivative_operator(1, p, xmin, xmax, N+1)
#         D2 = D^2
        D2 = periodic_derivative_operator(2, p, xmin, xmax, N+1)
        error_u = errors_bbm_periodic_manufactured(usol, D, D2, tspan, Tsit5(), tol, 1/N, true)
        push!(val_error, error_u)
    end
    Dict("val_p$(p)_N" => val_N, "val_p$(p)_error" => val_error)
end


filename = "bbm_periodic_manufactured_convergence_FD_narrow"
even_odd_values = evenodd_values
res = Dict{String,Any}()

p = 2
val_N = round.(Int, 2 .^ range(4, 10, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 4
val_N = round.(Int, 2 .^ range(4, 9, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 6
val_N = round.(Int, 2 .^ range(4, 8, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 8
val_N = round.(Int, 2 .^ range(4, 6.5, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

open("../data/" * filename * ".json", "w") do io
    JSON.print(io, res, 2)
end
```


```julia
fig, ax = plt.subplots(1, 1)
ax.set_prop_cycle(marker_cycler)

filename = "bbm_periodic_manufactured_convergence_FD_narrow"
res = JSON.parsefile("../data/" * filename * ".json")

let p = 2
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("\$ p = %d \$, EOC %.2f", p, -c1), (val_N[4], 2.0*val_error_u[4]), color="gray")
end

let p = 4
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("\$ p = %d \$, EOC %.2f", p, -c1), (val_N[6], 2.0*val_error_u[6]), color="gray")
end

let p = 6
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("\$ p = %d \$, EOC %.2f", p, -c1), (val_N[12], 2.0*val_error_u[12]), color="gray")
end

let p = 8
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("\$ p = %d \$, EOC %.2f", p, -c1), (0.15*val_N[end], val_error_u[end]), color="gray")
end

ax.set_ylim(1.0e-13, 1.5e0)
ax.set_xscale("log", basex=2)
ax.set_yscale("symlog", linthreshy=1.0e-13)
ax.set_xlabel(L"N")
ax.set_ylabel(L"\| u - u_{\mathrm{ana}} \|")
ax.locator_params(axis="y", numticks=5)
fig.savefig("../figures/" * filename * ".pdf", bbox_inches="tight")
```


```julia
function do_stuff(p, val_N)
    val_error = Float64[]
    for N in val_N
        D = periodic_derivative_operator(1, p, xmin, xmax, N+1)
        D2 = D^2
#         D2 = periodic_derivative_operator(2, p, xmin, xmax, N+1)
        error_u = errors_bbm_periodic_manufactured(usol, D, D2, tspan, Tsit5(), tol, 1/N, true)
        push!(val_error, error_u)
    end
    Dict("val_p$(p)_N" => val_N, "val_p$(p)_error" => val_error)
end


filename = "bbm_periodic_manufactured_convergence_FD_wide"
even_odd_values = evenodd_values
res = Dict{String,Any}()

p = 2
val_N = round.(Int, 2 .^ range(4, 10, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 4
val_N = round.(Int, 2 .^ range(4, 9, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 6
val_N = round.(Int, 2 .^ range(4, 8, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 8
val_N = round.(Int, 2 .^ range(4, 6.5, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

open("../data/" * filename * ".json", "w") do io
    JSON.print(io, res, 2)
end
```


```julia
fig, ax = plt.subplots(1, 1)
ax.set_prop_cycle(marker_cycler)

filename = "bbm_periodic_manufactured_convergence_FD_wide"
res = JSON.parsefile("../data/" * filename * ".json")

let p = 2
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("\$ p = %d \$, EOC %.2f", p, -c1), (val_N[4], 2.0*val_error_u[4]), color="gray")
end

let p = 4
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("\$ p = %d \$, EOC %.2f", p, -c1), (val_N[6], 2.0*val_error_u[6]), color="gray")
end

let p = 6
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("\$ p = %d \$, EOC %.2f", p, -c1), (val_N[12], 2.0*val_error_u[12]), color="gray")
end

let p = 8
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("\$ p = %d \$, EOC %.2f", p, -c1), (0.15*val_N[end], val_error_u[end]), color="gray")
end

ax.set_ylim(1.0e-13, 1.5e0)
ax.set_xscale("log", basex=2)
ax.set_yscale("symlog", linthreshy=1.0e-13)
ax.set_xlabel(L"N")
ax.set_ylabel(L"\| u - u_{\mathrm{ana}} \|")
ax.locator_params(axis="y", numticks=5)
fig.savefig("../figures/" * filename * ".pdf", bbox_inches="tight")
```

## Convergence study: Plots CG


```julia
function do_stuff(p, val_N)
    val_error = Float64[]
    for N in val_N
        mesh = UniformPeriodicMesh1D(xmin, xmax, N)
        D = couple_continuosly(legendre_derivative_operator(-1., 1., p+1), mesh)
        D2 = sparse(D)^2
#         D2 = sparse(couple_continuosly(legendre_second_derivative_operator(-1., 1., p+1), mesh))
        error_u = errors_bbm_periodic_manufactured(usol, D, D2, tspan, Tsit5(), tol, 1/N, true)
        push!(val_error, error_u)
    end
    Dict("val_p$(p)_N" => val_N, "val_p$(p)_error" => val_error)
end


filename = "bbm_periodic_manufactured_convergence_CG_D1c2"
even_odd_values = evenodd_values
res = Dict{String,Any}()

p = 1
val_N = round.(Int, 2 .^ range(6, 12, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 2
val_N = round.(Int, 2 .^ range(5, 10, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 3
val_N = round.(Int, 2 .^ range(4, 8, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 4
val_N = round.(Int, 2 .^ range(4, 7, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 5
val_N = round.(Int, 2 .^ range(2.5, 5.5, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 6
val_N = round.(Int, 2 .^ range(2, 4.75, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

open("../data/" * filename * ".json", "w") do io
    JSON.print(io, res, 2)
end
```


```julia
fig, ax = plt.subplots(1, 1)
ax.set_prop_cycle(marker_cycler)

filename = "bbm_periodic_manufactured_convergence_CG_D1c2"
res = JSON.parsefile("../data/" * filename * ".json")

let p = 1
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (val_N[2], 2.0*val_error_u[2]), color="gray")
end

let p = 2
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (val_N[end], 0.1*val_error_u[end]), color="gray")
end

let p = 3
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (val_N[end], 1.5*val_error_u[end]), color="gray")
end

let p = 4
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (0.4*val_N[end], 0.1*val_error_u[end]), color="gray")
end

let p = 5
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (1.2*val_N[end], 0.2*val_error_u[end]), color="gray")
end

let p = 6
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (0.2*val_N[end], 0.1*val_error_u[end]), color="gray")
end

plt.figure()
handles, labels = ax.get_legend_handles_labels()
plt.figlegend(handles, labels, loc="center", ncol=8)
plt.savefig("../figures/Galerkin_legend_p1_p6.pdf", bbox_inches="tight")

ax.set_xlim(2^1-0.5, 2^13-1)
ax.set_ylim(1.0e-13, 1.5e0)
ax.set_xscale("log", basex=2)
ax.set_yscale("symlog", linthreshy=1.0e-13)
ax.set_xlabel(L"N")
ax.set_ylabel(L"\| u - u_{\mathrm{ana}} \|")
ax.locator_params(axis="y", numticks=5)
fig.savefig("../figures/" * filename * ".pdf", bbox_inches="tight")
```


```julia
function do_stuff(p, val_N)
    val_error = Float64[]
    for N in val_N
        mesh = UniformPeriodicMesh1D(xmin, xmax, N)
        D = couple_continuosly(legendre_derivative_operator(-1., 1., p+1), mesh)
#         D2 = sparse(D)^2
        D2 = sparse(couple_continuosly(legendre_second_derivative_operator(-1., 1., p+1), mesh))
        error_u = errors_bbm_periodic_manufactured(usol, D, D2, tspan, Tsit5(), tol, 1/N, true)
        push!(val_error, error_u)
    end
    Dict("val_p$(p)_N" => val_N, "val_p$(p)_error" => val_error)
end


filename = "bbm_periodic_manufactured_convergence_CG_narrow"
even_odd_values = evenodd_values
res = Dict{String,Any}()

p = 1
val_N = round.(Int, 2 .^ range(6, 12, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 2
val_N = round.(Int, 2 .^ range(2, 8.5, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 3
val_N = round.(Int, 2 .^ range(2, 7, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 4
val_N = round.(Int, 2 .^ range(2, 5.5, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 5
val_N = round.(Int, 2 .^ range(2, 4.5, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 6
val_N = round.(Int, 2 .^ range(1, 4.25, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

open("../data/" * filename * ".json", "w") do io
    JSON.print(io, res, 2)
end
```


```julia
fig, ax = plt.subplots(1, 1)
ax.set_prop_cycle(marker_cycler)

filename = "bbm_periodic_manufactured_convergence_CG_narrow"
res = JSON.parsefile("../data/" * filename * ".json")

let p = 1
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (val_N[2], 2.0*val_error_u[2]), color="gray")
end

let p = 2
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (val_N[end], 2.0*val_error_u[end]), color="gray")
end

let p = 3
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (1.2*val_N[end], 0.8*val_error_u[end]), color="gray")
end

let p = 4
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (1.2*val_N[end], 0.2*val_error_u[end]), color="gray")
end

let p = 5
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (1.1*val_N[end], 0.1*val_error_u[end]), color="gray")
end

let p = 6
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (0.15*val_N[end], 0.5*val_error_u[end]), color="gray")
end

ax.set_xlim(2^1-0.5, 2^13-1)
ax.set_ylim(1.0e-13, 1.5e0)
ax.set_xscale("log", basex=2)
ax.set_yscale("symlog", linthreshy=1.0e-13)
ax.set_xlabel(L"N")
ax.set_ylabel(L"\| u - u_{\mathrm{ana}} \|")
ax.locator_params(axis="y", numticks=5)
fig.savefig("../figures/" * filename * ".pdf", bbox_inches="tight")
```

## Convergence study: Plots DG


```julia
function do_stuff(p, val_N)
    val_error = Float64[]
    for N in val_N
        Dop = legendre_derivative_operator(-1., 1., p+1); mesh = UniformPeriodicMesh1D(xmin, xmax, N)
        D = couple_discontinuosly(Dop, mesh); D₊ = couple_discontinuosly(Dop, mesh, Val(:plus)); D₋ = couple_discontinuosly(Dop, mesh, Val(:minus))
        D2 = sparse(D)^2
#         D2 = sparse(D₊) * sparse(D₋)
#         D2 = sparse(D₋) * sparse(D₊)
        error_u = errors_bbm_periodic_manufactured(usol, D, D2, tspan, Tsit5(), tol, 1/N, true)
        push!(val_error, error_u)
    end
    Dict("val_p$(p)_N" => val_N, "val_p$(p)_error" => val_error)
end


filename = "bbm_periodic_manufactured_convergence_DG_D1c2"
even_odd_values = evenodd_values
res = Dict{String,Any}()

p = 1
val_N = round.(Int, 2 .^ range(6, 12, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 2
val_N = round.(Int, 2 .^ range(5, 10, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 3
val_N = round.(Int, 2 .^ range(4, 8, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 4
val_N = round.(Int, 2 .^ range(3, 7, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 5
val_N = round.(Int, 2 .^ range(3, 5.75, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 6
val_N = round.(Int, 2 .^ range(2.5, 4.75, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

open("../data/" * filename * ".json", "w") do io
    JSON.print(io, res, 2)
end
```


```julia
fig, ax = plt.subplots(1, 1)
ax.set_prop_cycle(marker_cycler)

filename = "bbm_periodic_manufactured_convergence_DG_D1c2"
res = JSON.parsefile("../data/" * filename * ".json")

let p = 1
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (val_N[2], 2.0*val_error_u[2]), color="gray")
end

let p = 2
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (val_N[end], 2.0*val_error_u[end]), color="gray")
end

let p = 3
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (0.5*val_N[end], 0.2*val_error_u[end]), color="gray")
end

let p = 4
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (1.2*val_N[end], 0.2*val_error_u[end]), color="gray")
end

let p = 5
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (0.45*val_N[end], 0.1*val_error_u[end]), color="gray")
end

let p = 6
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (0.5*val_N[end], 0.1*val_error_u[end]), color="gray")
end

ax.set_xlim(2^2, 2^13-1)
ax.set_ylim(1.0e-13, 1.5e0)
ax.set_xscale("log", basex=2)
ax.set_yscale("symlog", linthreshy=1.0e-13)
ax.set_xlabel(L"N")
ax.set_ylabel(L"\| u - u_{\mathrm{ana}} \|")
ax.locator_params(axis="y", numticks=5)
fig.savefig("../figures/" * filename * ".pdf", bbox_inches="tight")
```


```julia
function do_stuff(p, val_N)
    val_error = Float64[]
    for N in val_N
        Dop = legendre_derivative_operator(-1., 1., p+1); mesh = UniformPeriodicMesh1D(xmin, xmax, N)
        D = couple_discontinuosly(Dop, mesh); D₊ = couple_discontinuosly(Dop, mesh, Val(:plus)); D₋ = couple_discontinuosly(Dop, mesh, Val(:minus))
#         D2 = sparse(D)^2
        D2 = sparse(D₊) * sparse(D₋)
#         D2 = sparse(D₋) * sparse(D₊)
        error_u = errors_bbm_periodic_manufactured(usol, D, D2, tspan, Tsit5(), tol, 1/N, true)
        push!(val_error, error_u)
    end
    Dict("val_p$(p)_N" => val_N, "val_p$(p)_error" => val_error)
end


filename = "bbm_periodic_manufactured_convergence_DG_D1pD1m"
even_odd_values = evenodd_values
res = Dict{String,Any}()

p = 1
val_N = round.(Int, 2 .^ range(6, 12, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 2
val_N = round.(Int, 2 .^ range(5, 10, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 3
val_N = round.(Int, 2 .^ range(4, 8, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 4
val_N = round.(Int, 2 .^ range(3, 7, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 5
val_N = round.(Int, 2 .^ range(3, 5.5, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 6
val_N = round.(Int, 2 .^ range(2.5, 4.5, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

open("../data/" * filename * ".json", "w") do io
    JSON.print(io, res, 2)
end
```


```julia
fig, ax = plt.subplots(1, 1)
ax.set_prop_cycle(marker_cycler)

filename = "bbm_periodic_manufactured_convergence_DG_D1pD1m"
res = JSON.parsefile("../data/" * filename * ".json")

let p = 1
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (val_N[2], 2.0*val_error_u[2]), color="gray")
end

let p = 2
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (val_N[end], 2.0*val_error_u[end]), color="gray")
end

let p = 3
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (0.5*val_N[end], 0.2*val_error_u[end]), color="gray")
end

let p = 4
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (1.2*val_N[end], 0.2*val_error_u[end]), color="gray")
end

let p = 5
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (0.7*val_N[end], 0.1*val_error_u[end]), color="gray")
end

let p = 6
    val_N, val_error_u = res["val_p$(p)_N"], res["val_p$(p)_error"]
    ax.plot(val_N, val_error_u, label="\$ p = $p \$")
    c0, c1 = linear_regression(log.(val_N), log.(val_error_u))
    plt.plot(val_N, exp(c0).*val_N.^c1, marker="", linestyle=":", color="gray")
    plt.annotate(@sprintf("EOC %.2f", -c1), (0.25*val_N[end], 0.1*val_error_u[end]), color="gray")
end

ax.set_xlim(2^2, 2^13-1)
ax.set_ylim(1.0e-13, 1.5e0)
ax.set_xscale("log", basex=2)
ax.set_yscale("symlog", linthreshy=1.0e-13)
ax.set_xlabel(L"N")
ax.set_ylabel(L"\| u - u_{\mathrm{ana}} \|")
ax.locator_params(axis="y", numticks=5)
fig.savefig("../figures/" * filename * ".pdf", bbox_inches="tight")
```


```julia

```


```julia

```

bd0f800653247aea7e542c39a5cc7797b63d3ef1
ipynb
Jupyter Notebook
notebooks/bbm.ipynb
ranocha/Dispersive-wave-schemes-notebooks
ff911d1451452f6667dc232597ce8870db28f627
MIT
notebooks/bbm.ipynb
ranocha/Dispersive-wave-schemes-notebooks
ff911d1451452f6667dc232597ce8870db28f627
MIT
notebooks/bbm.ipynb
ranocha/Dispersive-wave-schemes-notebooks
ff911d1451452f6667dc232597ce8870db28f627
MIT
Qwen/Qwen-72B
1. YES
2. YES