# Help to formulate an ODE problem to use OrdinaryDiffEq

**URL:** <https://discourse.julialang.org/t/help-to-formulate-an-ode-problem-to-use-ordinarydiffeq/112494>\
**Category:** Modelling & Simulations\
**Tags:** ode, ordinarydiffeq\
**Created:** [April 3, 2024, 11:10pm UTC](https://discourse.julialang.org/t/help-to-formulate-an-ode-problem-to-use-ordinarydiffeq/112494 "2024-04-03T23:10:24Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Iddingsite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iddingsite/32/31840_2.png) [@Iddingsite](https://discourse.julialang.org/u/Iddingsite)\
**Post date:** [April 3, 2024, 11:10pm UTC](https://discourse.julialang.org/t/help-to-formulate-an-ode-problem-to-use-ordinarydiffeq/112494/1 "2024-04-03T23:10:24Z")

</div>

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”:

```julia
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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 4, 2024, 8:05am UTC](https://discourse.julialang.org/t/help-to-formulate-an-ode-problem-to-use-ordinarydiffeq/112494/2 "2024-04-04T08:05:04Z")

</div>

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)`.

> [@Iddingsite](#):
>
> What complicates that is that a(X,t)a(X,t)a(X,t) doesn’t have an analytical expression, so a(\frac{X^{n+1} + x^{n}}{2}, t^{n+\frac{1}{2}})a(Xn+1+xn2,tn+12)a(\frac{X^{n+1} + x^{n}}{2}, t^{n+\frac{1}{2}}) needs to be interpolated from a(X, t^n)a(X,tn)a(X, t^n) and a(X, t^{n+1})a(X,tn+1)a(X, t^{n+1}) in each iteration

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.

---

<div class="post-metadata">

**Author:** ![Iddingsite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iddingsite/32/31840_2.png) [@Iddingsite](https://discourse.julialang.org/u/Iddingsite)\
**Post date:** [April 4, 2024, 9:36pm UTC](https://discourse.julialang.org/t/help-to-formulate-an-ode-problem-to-use-ordinarydiffeq/112494/3 "2024-04-04T21:36:52Z")

</div>

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:

```julia
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()`.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 5, 2024, 2:56am UTC](https://discourse.julialang.org/t/help-to-formulate-an-ode-problem-to-use-ordinarydiffeq/112494/4 "2024-04-05T02:56:30Z")

</div>

Are you trying to do a Velocity Verlet?

---

<div class="post-metadata">

**Author:** ![Iddingsite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iddingsite/32/31840_2.png) [@Iddingsite](https://discourse.julialang.org/u/Iddingsite)\
**Post date:** [April 5, 2024, 8:31am UTC](https://discourse.julialang.org/t/help-to-formulate-an-ode-problem-to-use-ordinarydiffeq/112494/5 "2024-04-05T08:31:24Z")

</div>

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

 ![SL_sketch_test](https://global.discourse-cdn.com/julialang/original/3X/c/6/c61bbad6a586ac160b94797ad352acce920b5537.jpeg)