# Strange behavior of integrator interface DiffEq.jl

**URL:** <https://discourse.julialang.org/t/strange-behavior-of-integrator-interface-diffeq-jl/37146>\
**Category:** Numerics\
**Tags:** diffeq\
**Created:** [April 7, 2020, 7:11am UTC](https://discourse.julialang.org/t/strange-behavior-of-integrator-interface-diffeq-jl/37146 "2020-04-07T07:11:04Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [April 7, 2020, 7:11am UTC](https://discourse.julialang.org/t/strange-behavior-of-integrator-interface-diffeq-jl/37146/1 "2020-04-07T07:11:04Z")

</div>

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.

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

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

---

<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 7, 2020, 1:32pm UTC](https://discourse.julialang.org/t/strange-behavior-of-integrator-interface-diffeq-jl/37146/2 "2020-04-07T13:32:50Z")

</div>

> [@mleprovost](#):
>
> Code if I integrate multiple copies of the initial condition `u0` stacked in an array `ens` , it doesn’t give the same result:

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.

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [April 7, 2020, 3:43pm UTC](https://discourse.julialang.org/t/strange-behavior-of-integrator-interface-diffeq-jl/37146/3 "2020-04-07T15:43:58Z")

</div>

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

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

---

<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 7, 2020, 5:02pm UTC](https://discourse.julialang.org/t/strange-behavior-of-integrator-interface-diffeq-jl/37146/4 "2020-04-07T17:02:27Z")

</div>

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