Hi everyone,

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

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.