How to add extra variables when defining ODEProblem in DifferentialEquations?

New to Julia
,
1

I’m trying to define a system of differential equations using matrix notation. The issue is that this matrix varies depending on the problem and in all the examples that I’ve seen the function f that goes in the ODEProblem() function, only has the three variables f(u,p,t).

In the following example, I want to define a function for which I can change the value that defines the size of the identity matrix (A):

using DifferentialEquations  
using LinearAlgebrape or paste code here

function f(u,p,t)
  A = I(4)
  p[1].*(A*u).*(1/t^(1+p[2]))
end

u0    = [0.1,0.3,0.5,0.8]
p     = [5.5, 0.5];

As you can deduct, I’m using the matrix A to define 4 different state variables. Each state variable will be defined by the same parameters, but the initial state of each variable is different.

tspan = (3.0,30.0)
prob = ODEProblem(f,u0,tspan,p)

Age1 = collect(5:5:30)
sol = solve(prob, saveat=Age1); #you can plot this for clarity

Ideally, I want something like this:

n = 4
function f(u,p,t,n)
  A = I(n)
  p[1].*(A*u).*(1/t^(1+p[2]))
end

I know this doesn’t work, but any ideas of how I can modify the code so it admits different sizes for the matrix A?

Thanks!

2

Could you just make p[3] = n ?

3

You could create a function that returns a function:

julia> function generate_f(n)
           A = I(n)
           f(u, p, t) = p[1].*(A*u).*(1/t^(1+p[2]))
           return f
       end
generate_f (generic function with 1 method)

julia> f4 = generate_f(4)
(::var"#f#8"{Diagonal{Bool, Vector{Bool}}}) (generic function with 1 method)

julia> f4(inits[1:4], inits[(nplot+1):(nplot+2)], 5)
4-element Vector{Float64}:
 2.726706751006412
 2.7139936889300844
 2.742306107133945
 2.6816291968428603

Another way is to create a functor.

julia> struct MyFunctor
           A::Diagonal{Bool, Vector{Bool}}
           MyFunctor(n::Int) = new(I(n))
       end

julia> (mf::MyFunctor)(u, p, t) = p[1].*(mf.A*u).*(1/t^(1+p[2]))

julia> f4(inits[1:4], inits[(nplot+1):(nplot+2)], 5)
4-element Vector{Float64}:
 2.7531945173074637
 2.74035795804277
 2.768945408652123
 2.7076790709064475

A third way would be to use a constant global Ref:

julia> const A_ref = Ref(I(4))
Base.RefValue{Diagonal{Bool, Vector{Bool}}}(Bool[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1])

julia> function f(u, p, t)
           p[1].*(A_ref[]*u).*(1/t^(1+p[2]))
       end
f (generic function with 1 method)

julia> f(inits[1:4], inits[(nplot+1):(nplot+2)], 5)
4-element Vector{Float64}:
 2.7531945173074637
 2.74035795804277
 2.768945408652123
 2.7076790709064475
4

It’s the identity matrix, so why not replace A*u with u and save the allocation?

2 Likes
5

You’re right. I thought I needed the matrix so it will recognize that it was a system of equations. But just by defining u as a vector, it recognizes it as a system.
So, this is just the same, without needing to define the matrix:

function f(u,p,t)
  p[1].*(u).*(1/t^(1+p[2]))
end

Tank you :slight_smile:

1 Like
6

I found this question to be relevant again. Now for an external variable that changes with time for each one of the variables in the system. This is an example:

Let’s say I have this, which works fine:

using DifferentialEquations  
using LinearAlgebra

function f(u,p,t)
    (p[1]+p[2]).*(u).*(1/t^(1+p[3]))
  end

u0    = [0.1,0.3,0.5,0.8]
p     = [1,3, 0.5];

tspan   = (3.0,30.0)
Age     = collect(5:5:30)
ndata   = length(Age)*length(u0)

prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob, saveat=Age)
plot(sol)

And I want to add the influence of another variable which varies for each one of the equations of the system and for each time for which we are solving the problem, something like this:

ext_var = rand(Float64, (length(u0),length(Age))) #matrix 

function g(u,p,t)
    (p[1]+p[2].*ext_var).*(u).*(1/t^(1+p[3])) #added extra var
  end

probg = ODEProblem(g,u0,tspan,p)
solg = solve(probg, saveat=Age)
plot(solg)

How should I define this problem. What I thought could work didn’t work.

7

All external parameters and data can be squeezed into p. The parameter argument p doesn’t have to be a vector - try using a tuple of parameters (including ext_var) instead. Since your ext_var parameter varies in time, you’ll need to use a callback to modify it. Following this example, try something like this:

p = (1, 3, 0.5,
    collect(eachcol(ext_var)), # time-varying parameter
    Ref(1))  # Ref container for ext_var column index

function g(u,p,t)
    ext_vec = p[4][p[5][]]
    @. (p[1]+p[2]*ext_vec)*(u)*(1/t^(1+p[3]))
end

# at the start of each age, increment the column index
affect!(integrator) = integrator.p[5][] += 1
cb = PresetTimeCallback(Age,affect!, save_positions=(false, false))

probg = ODEProblem(g,u0,tspan,p)
solg = solve(probg, callback=cb, saveat=Age)

I’m using the messy Ref[] indexing because tuples are immutable, but they can contain mutable elements, and Ref is a good single-element mutable container.

8

As already suggested, for numerical parameters you can/should use p. While non-array types are technically possible, this will cause problems with basically any and all parameter estimation and sensitivity analysis tools, so I would suggest ComponentArrays.

If you find yourself wanting to include other kinds of configuration options or static information, use a callable struct:

struct MyODE{T}
    opts::T
end

function (f::MyODE)(u,p,t)
    # do something with f.opts (or use it for dispatch)
    # ...
end

Edit: and for more complex use-cases with different model configurations, you might be interested in ModelParameters

9

an “external variable that changes with time” sounds like a forcing function, for which there is a standard approach.