Modifying time in sol after solve() function using DifferentialEquations.jl

Specific Domains Modelling & Simulations
, , , ,
1

Hi,

I was wondering, using DifferentialEquations.jl (and co), if there is an easy way to modify the t in sol object after it is output and to synchronise it with the interpolation?

I have tried to simply modify it with sol.t .=* 2 for instance, but it doesn’t seem to impact sol(x) when I try to obtain the values by interpolation.

If it is not clear, here is a MWE:

using OrdinaryDiffEq

function lorenz(u, p, t)
    dx = 10.0 * (u[2] - u[1])
    dy = u[1] * (28.0 - u[3]) - u[2]
    dz = u[1] * u[2] - (8 / 3) * u[3]
    [dx, dy, dz]
end

u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz, u0, tspan)
sol = solve(prob, Tsit5());

sol2 = copy(sol)

sol2.t .*= 2 

sol2(200) == sol(100)

Which returns false. I would like it to be true for my purpose.

I know it would be probably a bit of a hack, but my goal is to rescale back the values of sol.t to a different unit for instance. I know that I could just access that value by converting, but for my purpose, I would like to hide from the user the scaling part and just provide a solution object with the units he defines as input in another function.

2

Note you’re hitting internals here so use at your own risk. You would need to modify the internal interpolation time, which is sol2.interp.ts .*= 2. This is an internal detail so there is no guarantee this will continue to work, but I also don’t have reason to change it any time soon.

1 Like
3

Thanks! Exactly what I was looking for…

4

Well, interestingly, it actually doesn’t work properly. And I don’t really understand why:

using OrdinaryDiffEq

function lorenz(u, p, t)
    dx = 10.0 * (u[2] - u[1])
    dy = u[1] * (28.0 - u[3]) - u[2]
    dz = u[1] * u[2] - (8 / 3) * u[3]
    [dx, dy, dz]
end

u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz, u0, tspan)
sol = solve(prob, Tsit5());

sol2 = deepcopy(sol)

sol2.interp.ts .*= 2

@show sol2.t[end] == sol.t[end] * 2
# True
@show sol2(sol2.t[end]) == sol(sol.t[end])
# False
@show sol2(sol2.t[end])
# 3-element Vector{Float64}:
#  -4.361553927469487
#  -1.2879782191167406
#  26.917854626761358
@show sol(sol.t[end])
# 3-element Vector{Float64}:
#  -5.248620965966837
#  -1.4596434517264187
#  28.39189720840952

So now, sol.t got updated correctly, as does sol.interp.ts, but the interpolation is wrong. I don’t know if it is a rounding error or whatever with the timestep, but the last point is a modelled one, so it should be the same in both case.
I guess I shouldn’t play with that!