# 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}}: