DiffEq-Specific Array Types

In many cases, a standard array may not be enough to fully hold the data for a model. Many of the solvers in DifferentialEquations.jl (only the native Julia methods) allow you to solve problems on AbstractArray types which allow you to extend the meaning of an array. This page describes some of the AbstractArray types which can be helpful for modeling differential equations problems.

ArrayPartitions

ArrayPartitions in DiffEq are used for heterogeneous arrays. For example, DynamicalODEProblem solvers use them internally to turn the separate parts into a single array. You can construct an ArrayPartition using RecursiveArrayTools.jl:

using RecursiveArrayTools
A = ArrayPartition(x::AbstractArray...)

where x is an array of arrays. Then, A will act like a single array, and its broadcast will be type stable, allowing for it to be used inside of the native Julia DiffEq solvers in an efficient way. This is a good way to generate an array which has different units for different parts, or different amounts of precision.

Usage

An ArrayPartition acts like a single array. A[i] indexes through the first array, then the second, etc. all linearly. But A.x is where the arrays are stored. Thus for

using RecursiveArrayTools
A = ArrayPartition(y,z)

We would have A.x[1]==y and A.x[2]==z. Broadcasting like f.(A) is efficient.

Example: Dynamics Equations

In this example we will show using heterogeneous units in dynamics equations. Our arrays will be:

using Unitful, RecursiveArrayTools, DiffEqBase, OrdinaryDiffEq
using LinearAlgebra

r0 = [1131.340, -2282.343, 6672.423]u"km"
v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
Δt = 86400.0*365u"s"
μ = 398600.4418u"km^3/s^2"
rv0 = ArrayPartition(r0,v0)

Here, r0 is the initial positions, and v0 are the initial velocities. rv0 is the ArrayPartition initial condition. We now write our update function in terms of the ArrayPartition:

function f(dy, y, μ, t)
    r = norm(y.x[1])
    dy.x[1] .= y.x[2]
    dy.x[2] .= -μ .* y.x[1] / r^3
end

Notice that y.x[1] is the r part of y, and y.x[2] is the v part of y. Using this kind of indexing is type stable, even though the array itself is heterogeneous. Note that one can also use things like 2y or y.+x and the broadcasting will be efficient.

Now to solve our equations, we do the same thing as always in DiffEq:

prob = ODEProblem(f, rv0, (0.0u"s", Δt), μ)
sol = solve(prob, Vern8())

MultiScaleArrays

The multi-scale modeling functionality is provided by MultiScaleArrays.jl. It allows for designing a multi-scale model as an extension of an array, which in turn can be directly used in the native Julia solvers of DifferentialEquations.jl.

For more information, please see the MultiScaleArrays.jl README.

DEDataArrays

The DEDataArray{T} type allows one to add other "non-continuous" variables to an array, which can be useful in many modeling situations involving lots of events. To define an DEDataArray, make a type which subtypes DEDataArray{T} with a field x for the "array of continuous variables" for which you would like the differential equation to treat directly. The other fields are treated as "discrete variables". For example:

mutable struct MyDataArray{T,1} <: DEDataArray{T,1}
    x::Array{T,1}
    a::T
    b::Symbol
end

In this example, our resultant array is a SimType, and its data which is presented to the differential equation solver will be the array x. Any array which the differential equation solver can use is allowed to be made as the field x, including other DEDataArrays. Other than that, you can add whatever fields you please, and let them be whatever type you please.

These extra fields are carried along in the differential equation solver that the user can use in their f equation and modify via callbacks. For example, inside of a an update function, it is safe to do:

function f(du,u,p,t)
  u.a = t
end

to update the discrete variables (unless the algorithm notes that it does not step to the endpoint, in which case a callback must be used to update appropriately.)

Note that the aliases DEDataVector and DEDataMatrix cover the one and two dimensional cases.

Example: A Control Problem

In this example we will use a DEDataArray to solve a problem where control parameters change at various timepoints. First we will build

mutable struct SimType{T} <: DEDataVector{T}
    x::Array{T,1}
    f1::T
end

as our DEDataVector. It has an extra field f1 which we will use as our control variable. Our ODE function will use this field as follows:

function f(du,u,p,t)
    du[1] = -0.5*u[1] + u.f1
    du[2] = -0.5*u[2]
end

Now we will setup our control mechanism. It will be a simple setup which uses set timepoints at which we will change f1. At t=5.0 we will want to increase the value of f1, and at t=8.0 we will want to decrease the value of f1. Using the DiscreteCallback interface, we code these conditions as follows:

const tstop1 = [5.]
const tstop2 = [8.]


function condition(u,t,integrator)
  t in tstop1
end

function condition2(u,t,integrator)
  t in tstop2
end

Now we have to apply an effect when these conditions are reached. When condition is hit (at t=5.0), we will increase f1 to 1.5. When condition2 is reached, we will decrease f1 to -1.5. This is done via the functions:

function affect!(integrator)
  for c in full_cache(integrator)
    c.f1 = 1.5
  end
end

function affect2!(integrator)
  for c in full_cache(integrator)
    c.f1 = -1.5
  end
end

Notice that we have to loop through the full_cache array (provided by the integrator interface) to ensure that all internal caches are also updated. With these functions we can build our callbacks:

save_positions = (true,true)

cb = DiscreteCallback(condition, affect!, save_positions=save_positions)

save_positions = (false,true)

cb2 = DiscreteCallback(condition2, affect2!, save_positions=save_positions)

cbs = CallbackSet(cb,cb2)

Now we define our initial condition. We will start at [10.0;10.0] with f1=0.0.

u0 = SimType([10.0;10.0], 0.0)
prob = ODEProblem(f,u0,(0.0,10.0))

Lastly we solve the problem. Note that we must pass tstop values of 5.0 and 8.0 to ensure the solver hits those timepoints exactly:

const tstop = [5.;8.]
sol = solve(prob,Tsit5(),callback = cbs, tstops=tstop)

It's clear from the plot how the controls affected the outcome.

Data Arrays vs ParameterizedFunctions

The reason for using a DEDataArray is because the solution will then save the control parameters. For example, we can see what the control parameter was at every timepoint by checking:

[sol[i].f1 for i in 1:length(sol)]

A similar solution can be achieved using a ParameterizedFunction. We could have instead created our function as:

function f(du,u,p,t)
    du[1] = -0.5*u[1] + p
    du[2] = -0.5*u[2]
end
u0 = SimType([10.0;10.0], 0.0)
p = 0.0
prob = ODEProblem(f,u0,(0.0,10.0),p)
const tstop = [5.;8.]
sol = solve(prob,Tsit5(),callback = cbs, tstops=tstop)

where we now change the callbacks to changing the parameter:

function affect!(integrator)
  integrator.p = 1.5
end

function affect2!(integrator)
  integrator.p = -1.5
end

This will also solve the equation and get a similar result. It will also be slightly faster in some cases. However, if the equation is solved in this manner, there will be no record of what the parameter was at each timepoint. That is the tradeoff between DEDataArrays and ParameterizedFunctions.