# Fornberg-Whitham equation


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

# Periodic boundary conditions


```julia
function relaxation!(integrator)
    told = integrator.tprev
    uold = integrator.uprev
    tnew = integrator.t
    unew = integrator.u
    
    @unpack D, tmp1 = integrator.p
    
    @. tmp1 = (unew - uold)^2
    c = integrate(tmp1, D)
    @. tmp1 = (unew - uold) * uold
    γ = one(c)
    if !iszero(c)
        γ = -2 * integrate(tmp1, D) / c
    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
    
    nothing
end


function save_func_fw_periodic(u, t, integrator)
    @unpack D, invImD2, tmp1, usol = integrator.p
    
    mass = integrate(u, D)
    
    ldiv!(tmp1, invImD2, u)
    linear = integrate(tmp1, D)
    
    @. tmp1 = u^2
    quadratic = integrate(tmp1, D)
    
    x = grid(D)
    tmp1 .= ( usol.(t, x, x[1], -x[1]) .- u ).^2
    error_u = integrate(tmp1, D) |> sqrt
    
    SVector(mass, linear, quadratic)
end


function fw_periodic!(du, u, param, t)
    @unpack D, invImD2, tmp1, tmp2 = param
    one_third = one(t) / 3
    
    # energy and mass conservative split form semidiscretization
    mul!(tmp1, D, u)
    ldiv!(du, invImD2, tmp1)
    @. du += one_third * (u * tmp1)
    @. tmp1 = u * u
    mul!(tmp2, D, tmp1)
    @. du = -du - one_third * tmp2
    
    # standard central method without splitting: not energy conservative
#     mul!(tmp1, D, u)
#     ldiv!(du, invImD2, tmp1)
#     @. tmp1 = u * u
#     mul!(du, D, tmp1, -one(eltype(D))/2, -one(eltype(D)))
    
    nothing
end

function solve_ode_fw_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(fw_periodic!, u0, tspan, param)
    
    saveat = range(tspan..., length=100)
    saved_values_baseline = SavedValues(eltype(D), SVector{3,eltype(D)})
    saving_baseline = SavingCallback(save_func_fw_periodic, saved_values_baseline, saveat=saveat)
    saved_values_relaxation = SavedValues(eltype(D), SVector{3,eltype(D)})
    saving_relaxation = SavingCallback(save_func_fw_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)
    linear_baseline      = map(x->x[2], saved_values_baseline.saveval)
    linear_relaxation    = map(x->x[2], saved_values_relaxation.saveval)
    quadratic_baseline   = map(x->x[3], saved_values_baseline.saveval)
    quadratic_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)")
    plt.plot(t_relaxation, mass_relaxation .- mass_relaxation[1], label=L"$\int u$ (relaxation)")
    plt.plot(t_baseline,   linear_baseline   .- linear_baseline[1], label=L"$\int (u - \partial_x^2 u)$ (baseline)")
    plt.plot(t_relaxation, linear_relaxation .- linear_relaxation[1], label=L"$\int (u - \partial_x^2 u)$ (relaxation)")
    plt.plot(t_baseline,   quadratic_baseline    .- quadratic_baseline[1], label=L"$\int u^2$ (baseline)")
    plt.plot(t_relaxation, quadratic_relaxation  .- quadratic_relaxation[1], label=L"$\int u^2$ (relaxation)")
    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
# # single peakon solution
# xmin = -60.
# xmax = -xmin
# get_c() = 1.2
# function usol(t, x, xmin, xmax)
#     c = get_c()
#     x_t = mod(x - c*t - xmin, xmax - xmin) + xmin
    
#     4/3 * exp(-0.5 * abs(x_t)) + c - 4/3
# end

# traveling wave solution obtained numerically
c, data = open("fw_traveling_wave_init_c12_l80_N65536.txt", "r") do io
    line = readline(io)
    line = readline(io)
    line = readline(io)
    c = parse(Float64, line[8:end])
    data = readdlm(io, comments=true)
    c, data
end
x = range(data[1,1], data[end,1], length=size(data,1))
xmin = x[1]; xmax = x[end]
u0 = data[:, 2]
u0itp = CubicSplineInterpolation((x,), u0, extrapolation_bc=Periodic())
get_c() = c
function usol(t, x, xmin, xmax)
    c = get_c()
    x_t = mod(x - c*t - xmin, xmax - xmin) + xmin
    u0itp(x_t)
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())

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

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 = 8; 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 = 5; 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_fw_periodic(usol, D, D2, tspan, RK4(), tol, dt, adaptive)
```


```julia

```

# 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.diff(u^2/2, x, 3) |> 
        sp.simplify |> sp.string |> math_replacements |> println
end
```


```julia
function save_func_fw_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 fw_periodic_manufactured!(du, u, param, t)
    @unpack D, invImD2, tmp1, tmp2 = param
    x = grid(D)
    one_third = one(t) / 3
    
    mul!(tmp1, D, u)
    @. tmp2 = -tmp1 - 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)) - 16*pi^3*exp(t)*sinpi((2*t - 4*x)) - pi*exp(t)*sinpi((2*t - 4*x))
    ldiv!(du, invImD2, tmp2)
    @. du -= one_third * (u * tmp1)
    @. tmp1 = u * u
    mul!(tmp2, D, tmp1)
    @. du -= one_third * tmp2
    
    nothing
end

function errors_fw_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(fw_periodic_manufactured!, u0, tspan, param)
    
    saveat = range(tspan..., length=2)
    saved_values_baseline = SavedValues(eltype(D), eltype(D))
    saving_baseline = SavingCallback(save_func_fw_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-13

# 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 = 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_fw_periodic_manufactured(usol, D, D2, tspan, Tsit5(), tol, 1/N, true)
    push!(val_error, error_u)
end

fig, ax = plt.subplots(1, 1)
ax.set_prop_cycle(marker_cycler)
plt.plot(val_N, val_error, label=L"$\| u - u_{\mathrm{ana}} \|$")
c0, c1 = linear_regression(log.(val_N), log.(val_error))
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-12)
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_fw_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 = "fw_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, 10, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

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

p = 8
val_N = round.(Int, 2 .^ range(4, 7.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 = "fw_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), (0.6*val_N[end], 0.1*val_error_u[end]), 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}} \|")
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_fw_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 = "fw_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, 10, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

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

p = 8
val_N = round.(Int, 2 .^ range(4, 7.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 = "fw_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), (0.6*val_N[end], 0.1*val_error_u[end]), 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}} \|")
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_fw_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 = "fw_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(3, 6, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 6
val_N = round.(Int, 2 .^ range(3, 5.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 = "fw_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), (1.2*val_N[end], 0.5*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

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}} \|")
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_fw_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 = "fw_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(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(3, 6, length=7)) |> even_odd_values
res = merge(res, do_stuff(p, val_N))

p = 6
val_N = round.(Int, 2 .^ range(3, 5.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 = "fw_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], 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), (1.2*val_N[end], 0.5*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

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}} \|")
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_fw_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 = "fw_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(4, 7.15, 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(3, 5.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 = "fw_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.4*val_N[end], 0.1*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.6*val_N[1], 1.5*val_error_u[1]), 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}} \|")
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_fw_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 = "fw_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(4, 7.15, 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(3, 5.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 = "fw_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.4*val_N[end], 0.1*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.6*val_N[1], 1.5*val_error_u[1]), 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}} \|")
fig.savefig("../figures/" * filename * ".pdf", bbox_inches="tight")
```


```julia

```

6186b47407e62a293a55fc82d08c1f7c31968a04
ipynb
Jupyter Notebook
notebooks/fw.ipynb
ranocha/Dispersive-wave-schemes-notebooks
ff911d1451452f6667dc232597ce8870db28f627
MIT
notebooks/fw.ipynb
ranocha/Dispersive-wave-schemes-notebooks
ff911d1451452f6667dc232597ce8870db28f627
MIT
notebooks/fw.ipynb
ranocha/Dispersive-wave-schemes-notebooks
ff911d1451452f6667dc232597ce8870db28f627
MIT
Qwen/Qwen-72B
1. YES
2. YES