Help to formulate an ODE problem to use OrdinaryDiffEq

Hi everyone,

For one of my problem, I need to solve this ODE backwards:

\frac{dX}{dt} = a(X(t),t), \\ X(t^{n+1}) = X_g

where X is the position of a marker (could be 1d, 2d or 3d), a(X(t),t) is a velocity vector and X_g is a grid point.
The goal is to find the position of X at the previous timestep n knowing its current position at n+1 following characteristics and knowing a(X(t),t) at n and n+1 (this is for solving an advection problem using a semi-Lagrangian approach).

It is common to solve this equation using the implicit mid-point rule:

X(t^n) = - a(\frac{X(t^{n+1}) + X(t^n)}{2}, t^{n+\frac{1}{2}}) \Delta t + X(t^{n+1})

and iterating until convergence is reached as x^n is on both sides.
What complicates that is that a(X,t) doesn’t have an analytical expression, so a(\frac{X^{n+1} + x^{n}}{2}, t^{n+\frac{1}{2}}) needs to be interpolated from a(X, t^n) and a(X, t^{n+1}) in each iteration. So the final expression is

X(t^n)^{k+1} = - \mathcal{L}(\frac{X_g + X(t^n)^k}{2}, t^{n+\frac{1}{2}}) \Delta t + X_g, \: \textrm{for} \: k \ge 0, X(t^n)^{0} = X_g
with \mathcal{L} being a first order Lagrange interpolant for example.

So, I know how to program that, no problem. My issue is that I don’t know how I could combine that with OrdinaryDiffEq. My goal to do so is to try higher order method such as the implicit Adams–Moulton scheme and that would be nice to not have to implement that myself.

Here is a MWE of the approach “by hand”:

using Plots
using LinearAlgebra: norm

# grid
nx::Int = 51
Lx::Float64 = 1.0
Δx::Float64 = Lx / (nx-1)
x::StepRangeLen = range(0, length=nx, stop= Lx)
# velocity
v::Vector{Float64} = ones(nx) .* cos.(range(0, length=nx, stop= 10*Lx) .- 0.5) .+ 1.0
v_timestep::Vector{Float64} = ones(nx) .* cos.(range(0, length=nx, stop= 10*Lx) .- 0.2) .+ 1.0
v_t_half::Vector{Float64} = zeros(nx)
# CFL condition
CFL::Float64 = 1.5
period::Int = 2
tmax::Float64 = period * Lx / maximum(abs.(v))
Δt::Float64 = CFL * Δx / maximum(abs.(v))
X_t0::Vector{Float64} = zeros(nx)
X_t_half::Vector{Float64} = zeros(nx)

function linear_interpolation(x, y, new_x)
    idx = searchsortedlast(x, new_x)
    if idx == 0
        return y[1]
    elseif idx == length(x)
        return y[end]
    else
        frac = (new_x - x[idx]) / (x[idx+1] - x[idx])
        return y[idx] + frac * (y[idx+1] - y[idx])
    end
end

plt = scatter(x, zeros(size(x)))

plot!(plt, x, v, label="v")
plot!(plt, x, v_timestep, label="v_timestep")

v_t_half .= (v .+ v_timestep) / 2

# save the initial position
X_t_old = zeros(size(x))
# iterate 5 times to improve approximation
for i in 1:5
    X_t_old = copy(X_t_half)

    X_t_half .= x .- v_t_half .* Δt ./ 2

    v_t_half .= [linear_interpolation(x, v_t_half, xi) for xi in X_t_half]

    display(scatter!(X_t_half[1:10:end], zeros(size(X_t_half[1:10:end])), label=""))

    # compute error between X_t_half and X_t_old
    error = norm(X_t_half - X_t_old)
    println("Error: $error")
end

I fail to convert this code to something compatible with OrdinaryDiffEq.

I guess it should take something in that form:

using OrdinaryDiffEq

p = (v=v, v_timestep=v_timestep)

function system!(du, u, p, t)
    v, v_timestep = p
    du .= v
end

u0 = collect(x)
tspan = (0.0, -Δt)
prob = ODEProblem(system!, u0, tspan, dt = -Δt)

But I don’t know how I could tell it to interpolate v or if I should maybe provide this as a parameter.

1 Like

You just write down a as the function to solve, and make tspan be something that goes from a large time to small, like (1.0,0.0).

What do you mean by this? How is a defined? Is it given by a lookup table or interpolation or lookup table? If so, just put the lookup table calculation into the rhs.

Thx for the reply!

Ok, after thinking about it for a bit, a is given by an interpolation.

I think I maybe managed to do what I wanted:

using OrdinaryDiffEq

# grid
nx::Int = 51
Lx::Float64 = 1.0
Δx::Float64 = Lx / (nx-1)
x::StepRangeLen = range(0, length=nx, stop= Lx)
# velocity
v::Vector{Float64} = ones(nx) .* cos.(range(0, length=nx, stop= 10*Lx) .- 0.5) .+ 1.0
v_timestep::Vector{Float64} = ones(nx) .* cos.(range(0, length=nx, stop= 10*Lx) .- 0.2) .+ 1.0
v_t_half::Vector{Float64} = zeros(nx)
# CFL condition
CFL::Float64 = 1.5
period::Int = 2
tmax::Float64 = period * Lx / maximum(abs.(v))
Δt::Float64 = CFL * Δx / maximum(abs.(v))
X_t0::Vector{Float64} = zeros(nx)
X_t_half::Vector{Float64} = zeros(nx)

function linear_interpolation(x, y, new_x)
    idx = searchsortedlast(x, new_x)
    if idx == 0
        return y[1]
    elseif idx == length(x)
        return y[end]
    else
        frac = (new_x - x[idx]) / (x[idx+1] - x[idx])
        return y[idx] + frac * (y[idx+1] - y[idx])
    end
end

v_t_half .= (v .+ v_timestep) / 2

u0 = collect(x)
tspan = (Δt, 0)

p = (v_t_half, x)

function system!(du, u, p, t)
    v_t_half, x_final = p

    for i in axes(u, 1)
        du[i] = linear_interpolation(x_final, v_t_half, (u[i] .+ x_final[i]) / 2)
    end

end

prob = ODEProblem(system!, u0, tspan, p)

@btime sol = solve(prob, ImplicitMidpoint(), dt = -Δt)
# 179.603 μs (240 allocations: 86.91 KiB)
@btime sol = solve(prob, AB3(), dt = -Δt)
# 8.567 μs (55 allocations: 9.52 KiB)

Is that it? is that written in a general way? The results seem to be the same, but it seems to me that writting (u[i] .+ x_final[i]) / 2 already implies that I am using the ImplicitMidpoint().

Are you trying to do a Velocity Verlet?

No, this is part of a backward semi-Lagrangian, an advection algorithm often used in atmospheric or plasma modelling. Here is a sketch I did to illustrate what this step is supposed to do, finding the location at t_n of a marker X that is on the grid at t_{n+1} using a_{n+\frac{1}{2}}: