# Setup

This notebook has been created using Julia v1.5. Run the cells below to setup the project.


```julia
import Pkg
Pkg.activate(".")
Pkg.instantiate()
```


```julia
using LinearAlgebra

import Jacobi
using PolynomialBases
using SummationByPartsOperators

import PyPlot; plt = PyPlot
using LaTeXStrings, PyCall
using Latexify
import SymPy; sp = SymPy

ticker = pyimport("matplotlib.ticker")

# line cycler adapted to colourblind people
cycler = pyimport("cycler")
cycler.cycler
colours = ["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]
linestyles = ["-", "--", "-.", ":", "-", "--", "-."]
markers = ["4", "2", "3", "1", "+", "x", "."]
colourblind_cycler = (cycler.cycler(color=colours) + cycler.cycler(linestyle=linestyles))
plt.rc("axes", prop_cycle=colourblind_cycler)

plt.rc("text", usetex=true)
plt.rc("text.latex", preamble="\\usepackage{newpxtext}\\usepackage{newpxmath}\\usepackage{commath}\\usepackage{mathtools}")
plt.rc("font", family="serif", size=18.)
plt.rc("savefig", dpi=200)
plt.rc("legend", fontsize="medium", fancybox=true, framealpha=0.5)

isdir("../figures") || mkdir("../figures")
```

# Computation of the Butcher Coefficients


```julia
SummationByPartsOperators.mass_matrix(D::LegendreDerivativeOperator) = Diagonal(D.Δx * D.basis.weights)


function butcher_coefficients(M, D, x, tL, tR)
    Δ = tR'*x - tL'*x
    
    # grid oscillations
    kernel = nullspace(M \ D' * M)
    if size(kernel, 2) > 1
        throw(ArgumentError("The derivative operator D must be nullspace consistent!"))
    end
    osc = vec(kernel)
    
    # solve the linear system -> least norm solution
    A = qr(D, Val(true)) \ (I - (osc * osc' * M) / (osc'*M*osc))
    # and set the initial values to zero
    for j in 1:size(A,2)
        A[:,j] .-= tL' * A[:,j]
    end
    A ./= Δ
    
    A, diag(M) ./ Δ, (x .- tL' * x) ./ Δ
end

function butcher_coefficients(basis::NodalBasis)
    M = Diagonal(basis.weights)
    D = basis.D
    x = basis.nodes
    tL = basis.interp_left
    tR = basis.interp_right
    
    butcher_coefficients(M, D, x, tL, tR)
end

function butcher_coefficients(D_op::SummationByPartsOperators.AbstractDerivativeOperator)
    M = mass_matrix(D_op)
    D = Matrix(D_op)
    x = grid(D_op)
    tL = zero(x); tL[1] = 1
    tR = zero(x); tR[end] = 1
    
    butcher_coefficients(M, D, x, tL, tR)
end
```

# Lobatto Legendre Nodes: Lobatto IIIA Methods


```julia
D_op = legendre_derivative_operator(0., 1., 2)
A, b, c = butcher_coefficients(D_op)
display(A)
# the analytical values
[0 0; 1/2 1/2]
```


```julia
D_op = legendre_derivative_operator(0., 1., 3)
A, b, c = butcher_coefficients(D_op)
display(A)
# the analytical values
[0 0 0; 5/24 1/3 -1/24; 1/6 2/3 1/6]
```


```julia
D_op = legendre_derivative_operator(0., 1., 4)
A, b, c = butcher_coefficients(D_op)
display(A)
# the analytical values
sqrt5 = sqrt(5)
[0 0 0 0; 
    (11+sqrt5)/120 (25-sqrt5)/120 (25-13sqrt5)/120 (-1+sqrt5)/120;
    (11-sqrt5)/120 (25+13sqrt5)/120 (25+sqrt5)/120 (-1-sqrt5)/120;
    1/12 5/12 5/12 1/12
]
```


```julia

```

# Gauss/Radau Legendre Nodes: None of the Classical RK Schemes


```julia
basis = GaussLegendre(1)
A, b, c = butcher_coefficients(basis)
A
```


```julia
basis = GaussLegendre(2)
A, b, c = butcher_coefficients(basis)
display(A)
display(b)
display(c)

display(A - [-12+10*sqrt(15) -48+16*sqrt(15) -48+10*sqrt(15); 
        45+10*sqrt(15) 16*sqrt(15) -45+10*sqrt(15);
        48+10*sqrt(15) 48+16*sqrt(15) 12+10*sqrt(15)] / (72 * sqrt(15))) # A
display(b - [5, 8, 5] / 18) # b
display(c - [5 - sqrt(15), 5, 5 + sqrt(15)] / 10) # c
```


```julia
# folowwing the derivation using this example

let T = sp.symbols("T", real=true)
    basis = GaussLegendre(2, sp.Sym)
#     basis.D * 2/T / (sqrt(sp.Sym(15))/3T) |> latexify |> println
#     basis.nodes |> latexify |> println
#     basis.weights |> latexify |> println
    tL = sp.simplify(basis.interp_left)
#     tL * 6 |> println
    tR = sp.simplify(basis.interp_right)
#     tR * 6 |> println
    
    sqrt3 = sqrt(sp.Sym(3))
    τ = sp.simplify.((basis.nodes .+ 1) .* T / 2)
#     τ |> latexify |> println
#     τ = [T/2 - T/2sqrt3, T/2, T/2 + T/2sqrt3]
    D = sqrt(sp.Sym(15))/(3T) * [-3 4 -1; -1 0 1; 1 -4 3]
    M = Diagonal(T/18 * [5, 8, 5])
#     sp.simplify.(D * τ.^0)
#     k = 2; sp.simplify.(D * τ.^k - k * τ.^(k-1))
#     k = 2; sp.simplify(dot(basis.interp_left, τ.^k))
#     k = 2; sp.simplify(dot(basis.interp_right, τ.^k))
#     sp.simplify.(M * D + D' * M - (tR * tR' - tL * tL'))
#     M \ (D' * M) * 30T/sqrt(sp.Sym(15))
#     sp.simplify.((3 .* τ.^2 .- 1) * (-5 / (3*T^2/4 - 1))) # highest resolvable Legendre polynomial
    osc = [4, -5, 4]
#     osc' * M * D
#     osc' * M
    F = sp.Sym.(I - osc * osc' * M / dot(osc, M, osc))
    F * 18
    
    sqrt15 = sqrt(sp.Sym(15))
    a1, a2, b1, b2 = sp.symbols("a_1, a_2, b_1, b_2", real=true)
    basis_V0 = [a1, a2, -(4+sqrt15)*a1 + 2*(5+sqrt15)*a2/5]
    basis_V1 = [b1, b2, -b1 + 2*b2]
    J = 3T / (24*sqrt15) * [2 4*(sqrt15-3) -2; 5 4sqrt15 -5; 2 4*(sqrt15+3) -2]
    sp.simplify.(tL' * J)
    sp.simplify.(J * D * 6)
    sp.simplify.(J * D * basis_V0)
    sp.simplify(D * J * 2)
    sp.simplify.(D * J * basis_V1)
    
    sp.simplify.(J * F / (T/72sqrt15))
    M
    τ * 10/T
end
```


```julia
q = Jacobi.Quadrature(Jacobi.GRJM, 2)
M = Diagonal(Jacobi.qweights(q))
D = Jacobi.qdiff(q)
x = Jacobi.qzeros(q)
tL = Jacobi.interp_mat([-one(eltype(x))], x) |> vec
tR = Jacobi.interp_mat([+one(eltype(x))], x) |> vec
A, b, c = butcher_coefficients(M, D, x, tL, tR)
A
```


```julia

```

# The FD SBP Operators used in the Numerical Tests


```julia
N = 3
order = 2
D_op = derivative_operator(MattssonNordström2004(), 1, order, 0//1, 1//1, N)
A_, b, c = butcher_coefficients(D_op)
A = rationalize.(A_, tol=1.0e-5)
display(A)
display(b)
display(c)
latexify(A) |> println
latexify(b) |> println
latexify(c) |> println
```


```julia
N = 9
order = 2
D_op = derivative_operator(MattssonNordström2004(), 1, order, 0//1, 1//1, N)
A_, b, c = butcher_coefficients(D_op)
A = rationalize.(A_, tol=1.0e-5)
display(A)
display(b)
display(c)
latexify(A) |> println
latexify(b) |> println
latexify(c) |> println
```


```julia
N = 9
order = 4
D_op = derivative_operator(MattssonNordström2004(), 1, order, 0//1, 1//1, N)
A_, b, c = butcher_coefficients(D_op)
A = rationalize.(A_, tol=1.0e-5)
display(A)
display(b)
display(c)
latexify(A) |> println
latexify(b) |> println
latexify(c) |> println
```


```julia

```


```julia

```

# Numerical Tests Using FD SBP Operators


```julia
function compute_error_nonstiff(λ, N, order)
    D_op = derivative_operator(MattssonNordström2004(), 1, order, 0., 1., N)
    
    M = mass_matrix(D_op)
    D = Matrix(D_op)
    t = SummationByPartsOperators.grid(D_op)
    Δ = last(t) - first(t)
    
    # grid oscillations
    kernel = nullspace(M \ D' * M)
    if size(kernel, 2) != 1
        throw(ArgumentError("The derivative operator D must be nullspace consistent!"))
    end
    osc = vec(kernel)
    
    # solve the linear system -> least norm solution
    A = qr(D, Val(true)) \ (I - (osc * osc' * M) / (osc'*M*osc))
    # and set the initial values to zero
    for j in 1:size(A,2)
        A[:,j] .-= A[1,j]
    end
    
    # compute solution
    u0 = fill(1., size(D,2))
    u = (I - λ*A) \ u0
    
    t, u, abs(last(u) - exp(λ*(last(t)-first(t))))
end

function compute_error_stiff(λ, N, order)
    D_op = derivative_operator(MattssonNordström2004(), 1, order, 0., 1., N)
    
    M = mass_matrix(D_op)
    D = Matrix(D_op)
    t = SummationByPartsOperators.grid(D_op)
    Δ = last(t) - first(t)
    
    # grid oscillations
    kernel = nullspace(M \ D' * M)
    if size(kernel, 2) != 1
        throw(ArgumentError("The derivative operator D must be nullspace consistent!"))
    end
    osc = vec(kernel)
    
    # solve the linear system -> least norm solution
    A = qr(D, Val(true)) \ (I - (osc * osc' * M) / (osc'*M*osc))
    # and set the initial values to zero
    for j in 1:size(A,2)
        A[:,j] .-= A[1,j]
    end
    
    # compute solution
    u0 = fill(1., size(D,2))
    u = (I - λ*A) \ (u0 - A * (@. (λ+1) * exp(-t)))
    
    t, u, abs(last(u) - exp(-(last(t)-first(t))))
end

function convergence_orders(Ns, errors)
    @assert length(Ns) == length(errors)
    
    orders = zero(errors)
    orders[1] = NaN
    for i in 2:length(errors)
        orders[i] = - log(errors[i-1] / errors[i]) / log(Ns[i-1] / Ns[i])
    end
    orders
end

#t, u, e = compute_error_stiff(-1000.0, 100, 2); println(e); Plots.plot(t, u); Plots.plot!(t, exp.(-t))
```

## Non-stiff Problem


```julia
Ns = 20:10:130

plt.close("all")
fig, ax = plt.subplots(1, 1)

for (idx,order) in enumerate(2:2:8)
    errors = zeros(length(Ns))
    for (idx,N) in enumerate(Ns)
        _, _, errors[idx] = compute_error_nonstiff(-1.0, N, order)
    end
    plt.plot(Ns, errors, label="Interior order $order", marker=markers[idx], ms=12, mew=1.5, lw=2)
    println(convergence_orders(Ns, errors))
    
    order_color = "silver"
    if order == 2
        plt.plot([Ns[1],Ns[end]], 7.e-5.*[1, (Ns[1]/Ns[end])^2], "-", color=order_color)
        plt.annotate("Order 2", (Ns[1], 6.e-6), color=order_color)
    elseif order == 4
        plt.plot([Ns[1],Ns[end]], 5.e-8.*[1, (Ns[1]/Ns[end])^4], "-", color=order_color)
        plt.annotate("Order 4", (Ns[1], 2.e-9), color=order_color)
    elseif order == 6
        plt.plot([Ns[1],Ns[end]], 1.e-10.*[1, (Ns[1]/Ns[end])^6], "-", color=order_color)
        plt.annotate("Order 6", (Ns[1], 1.5e-10), color=order_color)
    elseif order == 8
        plt.plot([Ns[1],60], 7.e-13.*[1, (Ns[1]/60)^8], "-", color=order_color)
        plt.annotate("Order 8", (Ns[1], 7.e-15), color=order_color)
    end
end
sleep(0.1)

plt.xlabel(L"N")
plt.xscale("log")
plt.ylabel("Error")
plt.yscale("log")
ax.xaxis.set_minor_formatter(ticker.NullFormatter())
ax.set_xticks([20, 30, 40, 60, 80, 100, 130], [])
ax.get_xaxis().set_major_formatter(ticker.ScalarFormatter())
plt.legend(loc="upper left", bbox_to_anchor=(1,1))
plt.savefig("../figures/convergence_scalar_nonstiff.pdf", bbox_inches="tight")
```

## Stiff Problem


```julia
Ns = 20:20:200

plt.close("all")
fig, ax = plt.subplots(1, 1)

for (idx,order) in enumerate(2:2:8)
    errors = zeros(length(Ns))
    for (idx,N) in enumerate(Ns)
        _, _, errors[idx] = compute_error_stiff(-1000.0, N, order)
    end
    plt.plot(Ns, errors, label="Interior order $order", marker=markers[idx], ms=12, mew=1.5, lw=2)
    println(convergence_orders(Ns, errors))
    
    order_color = "silver"
    if order == 2
        plt.plot([Ns[1],80], 2.e-5.*[1, (Ns[1]/80)^1], "-", color=order_color)
        plt.annotate("Order 1", (30, 1.5e-5), color=order_color)
        plt.plot([80,Ns[end]], 5.e-6.*[1, (80/Ns[end])^2], "-", color=order_color)
        plt.annotate("Order 2", (80, 6.e-6), color=order_color)
    elseif order == 4
        plt.plot([Ns[1],Ns[end]], 2.5e-6.*[1, (Ns[1]/Ns[end])^2.3], "-", color=order_color)
        plt.annotate("Order 2.3", (Ns[1], 3.e-6), color=order_color)
    elseif order == 6
        plt.plot([Ns[1],Ns[end]], 5.e-8.*[1, (Ns[1]/Ns[end])^3], "-", color=order_color)
        plt.annotate("Order 3", (Ns[1], 6.e-8), color=order_color)
    elseif order == 8
        plt.plot([Ns[1],Ns[end]], 7.e-9.*[1, (Ns[1]/Ns[end])^4.4], "-", color=order_color)
        plt.annotate("Order 4.4", (Ns[1], 2.5e-10), color=order_color)
    end
end
sleep(0.1)

plt.xlabel(L"N")
plt.xscale("log")
plt.ylabel("Error")
plt.yscale("log")
ax.xaxis.set_minor_formatter(ticker.NullFormatter())
ax.set_xticks([20, 30, 40, 60, 80, 100, 140, 200], [])
ax.get_xaxis().set_major_formatter(ticker.ScalarFormatter())
plt.legend(loc="upper left", bbox_to_anchor=(1,1))
plt.savefig("../figures/convergence_scalar_stiff.pdf", bbox_inches="tight")
```


```julia

```

eff2a9257bd2c19443c4f806d3194a4a4493f4fb
ipynb
Jupyter Notebook
notebooks/SBP_projection_time_integration.ipynb
ranocha/SBP-projection-in-time-notebooks
f97b25149698c8669afd158862cd76d35fbca181
MIT
notebooks/SBP_projection_time_integration.ipynb
ranocha/SBP-projection-in-time-notebooks
f97b25149698c8669afd158862cd76d35fbca181
MIT
notebooks/SBP_projection_time_integration.ipynb
ranocha/SBP-projection-in-time-notebooks
f97b25149698c8669afd158862cd76d35fbca181
MIT
Qwen/Qwen-72B
1. YES
2. YES