Strange behavior of integrator interface DiffEq.jl

Specific Domains Numerics
1

Hello,

I would like to use the integrator interface of DifferentialEquations.jl for an ensemble, i.e. loop over the different ensemble members and advance forward in time by one time step each of them.

However, if I integrate multiple copies of the initial condition u0 stacked in an array ens, it doesn’t give the same result.

using Test
using DifferentialEquations, OrdinaryDiffeq
function lorenz(du,u,p,t)
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end

u0 = [10.0; -5.0; 2.0]
tspan = (0.0,40.0)
Δt = 1e-2
T = tspan[1]:Δt:tspan[end]

prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob, RK4(), adaptive = false, dt = Δt);
 
# Integrator approach
integrator = init(prob, RK4(), adaptive =false, dt = Δt, save_everystep=false)

states = [deepcopy(u0)]
for t in T[1:end-1]
    step!(integrator)
    push!(states, deepcopy(integrator.u))
end

@test norm(states[end]-sol.u[end])<1e-15

Test Passed

Code if I integrate multiple copies of the initial condition u0 stacked in an array ens, it doesn’t give the same result:

k = 3
Ne = 5

ens = zeros(k,Ne)
ens .= repeat(u0, 1, Ne)
#Define propagation operator
function Propagate(t::Float64, ens, integrator::OrdinaryDiffEq.ODEIntegrator) 
    Nx, Ne = size(ens)
    for i=1:Ne
        s = deepcopy(ens[:,i])
        set_t!(integrator, deepcopy(t))
        set_u!(integrator, deepcopy(s))
        step!(integrator)
        ens[:,i] .= deepcopy(integrator.u)
    end
end

ens_hist = []
# Push the initial state
push!(ens_hist, deepcopy(ens))

# Propagate each ensemble member
for t in T[1:end-1]
    Propagate(t, ens, integrator)
    push!(ens_hist, deepcopy(ens))
end

@show sol.u[end]
@show ens_hist[end][:,1]

sol.u[end] = [-4.224313362730708, 0.6368037936621876, 28.848001043272014]

3-element Array{Float64,1}:
-4.224313362730708
0.6368037936621876
28.848001043272014

(ens_hist[end])[:, 1] = [-4.229515264992412, 0.6012730192194168, 28.81877934048208]

3-element Array{Float64,1}:
-4.229515264992412
0.6012730192194168
28.81877934048208

2

I’m not sure why this is strange. You’re doing set_t! and set_u! instead of reinit!, so your starting timestepping PI values are not reset. This means you get slightly different stepping behavior. If you do uncertainty quantification (can link after GH-pages is done with its maintenance) you can see that the exponential divergence due to chaos kicks in, so even a difference of abstol=1e-6 vs abstol=1f-6 is enough to give an O(1) difference.

1 Like
3

I am now using reinit! but I still get a different result.

k = 3
Ne = 5
# ens = EnsembleState(3, 10)
ens = zeros(k,Ne)
ens .= repeat(u0, 1, Ne)
#Define propagation operator
function Propagate(t::Float64, ens, integrator::OrdinaryDiffEq.ODEIntegrator) 
    Nx, Ne = size(ens)
    for i=1:Ne
        s = ens[:,i]
        reinit!(integrator, deepcopy(s), t0 = deepcopy(t))
#         set_t!(integrator, deepcopy(t))
#         set_u!(integrator, deepcopy(s))
#         @show integrator
        # Choose the number of time steps to advance the system
#         for j=1:1
        step!(integrator)
#         @show integrator.u
        ens[:,i] .= deepcopy(integrator.u)
    end
end

ens_hist = []
push!(ens_hist, deepcopy(ens))
for t in T[1:end-1]
    Propagate(t, ens, integrator)
    push!(ens_hist, deepcopy(ens))
end

@show sol.u[end]
@show ens_hist[end][:,1]

sol.u[end] = [-4.224313362730708, 0.6368037936621876, 28.848001043272014]

3-element Array{Float64,1}:
-4.224313362730708
0.6368037936621876
28.848001043272014

(ens_hist[end])[:, 2] = [-4.224313362701003, 0.6368037936639874, 28.84800104322336]

3-element Array{Float64,1}:
-4.224313362701003
0.6368037936639874
28.84800104322336

1 Like
4

that’s highly likely to just be round off in the 1e-16 decimal place.

1 Like