Handle differential equations with size varying state vector

General Usage
1

Hello,

I am working on time marching of ODEs with size varying state vector. I would like to use the integrator interface of DifferentialEquations to do so. Here is a toy problem with the error. du is computed based on the initial u and not on modified u.

using DifferentialEquations

function RHS(du,u,p,t)
    N = size(u,1)
    for i = 1:N
    du[i] = -u[i]
    end
@show size(du)
end

u0 = [1.0; 1.0; 1.0]
states = [deepcopy(u0)]
tspan = (0.0,1.0)

Δt = 1e-1
T = tspan[1]:Δt:tspan[end]

prob = ODEProblem(RHS,u0,tspan)

integrator = init(prob, RK4(), adaptive =false, dt = Δt, save_everystep=false)

returns


t: 0.0
u: 3-element Array{Float64,1}:
 1.0
 1.0
 1.0

For loop for the time marching

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

returns

size(du) = (3,)
size(du) = (3,)
size(du) = (3,)

DimensionMismatch("array could not be broadcast to match destination")

DimensionMismatch("array could not be broadcast to match destination")

Stacktrace:
 [1] check_broadcast_shape at ./broadcast.jl:503 [inlined]
 [2] check_broadcast_axes at /Users/mathieu/.julia/packages/DiffEqBase/E16PL/src/diffeqfastbc.jl:19 [inlined]
 [3] check_broadcast_axes at ./broadcast.jl:509 [inlined]
 [4] check_broadcast_axes at ./broadcast.jl:510 [inlined]
 [5] instantiate at ./broadcast.jl:259 [inlined]
 [6] materialize! at ./broadcast.jl:801 [inlined]
 [7] macro expansion at /Users/mathieu/.julia/packages/DiffEqBase/E16PL/src/diffeqfastbc.jl:79 [inlined]
 [8] perform_step!(::OrdinaryDiffEq.ODEIntegrator{RK4,true,Array{Float64,1},Float64,DiffEqBase.NullParameters,Float64,Float64,Float64,Array{Array{Float64,1},1},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},RK4,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.RK4Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1}}},DiffEqBase.DEStats},ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.RK4Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64,Nothing}, ::OrdinaryDiffEq.RK4Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1}}, ::Bool) at /Users/mathieu/.julia/packages/OrdinaryDiffEq/IjEvf/src/perform_step/fixed_timestep_perform_step.jl:320
 [9] perform_step! at /Users/mathieu/.julia/packages/OrdinaryDiffEq/IjEvf/src/perform_step/fixed_timestep_perform_step.jl:309 [inlined]
 [10] step!(::OrdinaryDiffEq.ODEIntegrator{RK4,true,Array{Float64,1},Float64,DiffEqBase.NullParameters,Float64,Float64,Float64,Array{Array{Float64,1},1},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},RK4,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.RK4Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1}}},DiffEqBase.DEStats},ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.RK4Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64,Nothing}) at /Users/mathieu/.julia/packages/OrdinaryDiffEq/IjEvf/src/iterator_interface.jl:12
 [11] top-level scope at ./In[17]:3
2

Does this fit your problem:
https://docs.juliadiffeq.org/dev/features/callback_functions/#Example-3:-Growing-Cell-Population-1

?

4 Likes
3

Thanks for your answer @mauro3