How to add extra variables when defining ODEProblem in DifferentialEquations?

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!

Could you just make p[3] = n ?

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

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

2 Likes

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