How to define a system of ODEs with @ode_def and update the value of a variable using a function that depends on time?

New to Julia
, , ,
1

Hello, I’m new to Julia and I have a problem with my code. I am defining a system of differential equations using @ode_def of ParameterizedFunctions.jl and I want to determine the value of a variable using a function. However, I only want to change the value of the variable if t > 0. Doing this operation gives me an error because I am using the symbol t in a boolean context.

I know I could do this with a callback, but even though it doesn’t throw errors, it doesn’t give me the correct result.

What should I do? Is there a way to continue using @ode_def or should I simply use a generic function to define the system? I would prefer to find a solution to continue using @ode_def.

Here is the code:

using DifferentialEquations
using ParameterizedFunctions

# Solve odeFB
odeFB = @ode_def begin 
    dGp = ((Vmax * Ge * Gp)/(Km + Gp)) - (k4 * Gp) - (k2 * A * Gp) + (km2 * A) + (km4 * Gi)
    dGi = (k4 * Gp) - (k5 * Gi) - (km4 * Gi) + (km5 * A)
    dA = - (k3 * A) + (k5 * Gi) - (km5 * A)
end Vmax Ge Km k2 k3 k4 k5 km2 km4 km5

x0 = ones(length(odeFB.syms))
p  = (1400, 24.15, 1100, 5.3, 4, 2, 2.3, 5.3, 2, 2.3)
solFB = solve(ODEProblem(odeFB,x0,1000,p); reltol = 1e-12);

# External function to calculate a constant of the equation du[2] of the second system 
function f(t, solFB)
    if t > 0
        i = findfirst(times .>= t)
        return ((((solFB[1,i] - solFB[1,i-1]) / (solFB.t[i] - solFB.t[i-1])) * (t .- solFB.t[i-1])) .+ solFB[1,i-1])
    else 
        return solFB[1,1]
    end
end 

# Solve odeNF
odeNF = @ode_def begin 
    dGp = ((Vmax * Ge * Gp)/(Km + Gp)) - (k4 * Gp) - (k2 * A * Gp) + (km2 * A) + (km4 * Gi)
    dGi = (k4*f(t, solFB)) - (k5 * Gi) - (km4 * Gi) + (km5 * A)
    dA = - (k3 * A) + (k5 * Gi) - (km5 * A)
end Vmax Ge Km k2 k3 k4 k5 k6 km2 km4 km5 

solNF = solve(ODEProblem(odeNF, x0, 1000, p); reltol = 1e-12);
2

You need to @register_symbolic it so that way it’s kept as a node in the symbolic graph.

using DifferentialEquations
using ParameterizedFunctions

# Solve odeFB
odeFB = @ode_def begin 
    dGp = ((Vmax * Ge * Gp)/(Km + Gp)) - (k4 * Gp) - (k2 * A * Gp) + (km2 * A) + (km4 * Gi)
    dGi = (k4 * Gp) - (k5 * Gi) - (km4 * Gi) + (km5 * A)
    dA = - (k3 * A) + (k5 * Gi) - (km5 * A)
end Vmax Ge Km k2 k3 k4 k5 km2 km4 km5

x0 = ones(length(odeFB.syms))
p  = (1400, 24.15, 1100, 5.3, 4, 2, 2.3, 5.3, 2, 2.3)
solFB = solve(ODEProblem(odeFB,x0,1000,p); reltol = 1e-12);

# External function to calculate a constant of the equation du[2] of the second system 
function f(t, solFB)
    if t > 0
        i = findfirst(solFB.t .>= t)
        return ((((solFB[1,i] - solFB[1,i-1]) / (solFB.t[i] - solFB.t[i-1])) * (t .- solFB.t[i-1])) .+ solFB[1,i-1])
    else 
        return solFB[1,1]
    end
end
f2(t) = f(t, solFB)

using Symbolics
@register_symbolic f2(t)

# Solve odeNF
odeNF = @ode_def begin 
    dGp = ((Vmax * Ge * Gp)/(Km + Gp)) - (k4 * Gp) - (k2 * A * Gp) + (km2 * A) + (km4 * Gi)
    dGi = (k4 * f2(t)) - (k5 * Gi) - (km4 * Gi) + (km5 * A)
    dA = - (k3 * A) + (k5 * Gi) - (km5 * A)
end Vmax Ge Km k2 k3 k4 k5 km2 km4 km5
prob = ODEProblem(odeNF, x0, 1000, p); 
prob.f(x0,p,0.0)

solNF = solve(prob, reltol = 1e-12);

See Function Registration and Tracing · Symbolics.jl for more details.

This is because ParameterizedFunctions.jl is implemented using ModelingToolkit.jl/Symbolics.jl, and this is a property of the symbolic library. Generally I wouldn’t recommend ParameterizedFunctions.jl that much these days except for very simple cases, but you can use this to make it do crazy things just by registering more functions to it.

1 Like
3

It works perfectly! Thank you very much! Just out of curiosity, why wouldn’t you recommend using ParameterizedFunctions.jl anymore?

4

Ehh it’s just very limited in comparison to ModelingToolkit without too many gains. People hit the wall of what it can do and ask for more features all of the time, and ModelingToolkit is the more flexible version. Nothing wrong with ParameterizedFunctions if you’re keeping it to simple things like this though.

1 Like
5

Hey! I was playing with your solution, and it works perfectly until you change the second parameter from 24.15 to 18 (just to give an example; there is a whole range of parameters it does not work with). The error that I get is the following:

ERROR: LoadError: MethodError: Cannot `convert` an object of type SymbolicUtils.BasicSymbolic{Float64} to an object of type Float64

It’s really odd.

Anyway, thanks for helping me before!

6

That shouldn’t do anything, the exact numbers don’t matter in the tracing.

7

Yes, I’m quite sure it has to do with the ODE solver, as the error doesn’t occur with certain algorithm, such as KenCarp4().

8

This would require an MWE.

9

Okay, here’s a MWE where a system of ODEs uses an external function that will always return zero to add it to its first equation, therefore the external function shouldn’t affect the solution. The first example is identical to the second, with the exception that in the second example, it’s specified that KenCarp4() should be used to solve the system.

Throws an error 
using DifferentialEquations
using ParameterizedFunctions

# External function
function f(t)
    return t*0
end

using Symbolics
@register_symbolic f(t)

# ODE
ode = @ode_def begin 
    dGp = ((Vmax * Ge * Gp)/(Km + Gp)) - (k4 * Gp) - (k2 * A * Gp) + (km2 * A) + (km4 * Gi) + f(t)
    dGi = (k4 * Gp) - (k5 * Gi) - (km4 * Gi) + (km5 * A)
    dA = - (k3 * A) + (k5 * Gi) - (km5 * A)
end Vmax Ge Km k2 k3 k4 k5 km2 km4 km5

x0 = ones(length(ode.syms))
p  = (1400, 18, 1100, 5.3, 4, 2, 2.3, 5.3, 2, 2.3)
prob = ODEProblem(ode, x0, 1000, p); 

sol = solve(prob, reltol = 1e-12);

Error:

ERROR: LoadError: MethodError: Cannot `convert` an object of type SymbolicUtils.BasicSymbolic{Float64} to an object of type Float64
Works fine 
using DifferentialEquations
using ParameterizedFunctions

# External function 
function f(t)
    return t*0
end

using Symbolics
@register_symbolic f(t)

# ODE
ode = @ode_def begin 
    dGp = ((Vmax * Ge * Gp)/(Km + Gp)) - (k4 * Gp) - (k2 * A * Gp) + (km2 * A) + (km4 * Gi) + f(t)
    dGi = (k4 * Gp) - (k5 * Gi) - (km4 * Gi) + (km5 * A)
    dA = - (k3 * A) + (k5 * Gi) - (km5 * A)
end Vmax Ge Km k2 k3 k4 k5 km2 km4 km5

x0 = ones(length(ode.syms))
p  = (1400, 18, 1100, 5.3, 4, 2, 2.3, 5.3, 2, 2.3)
prob = ODEProblem(ode, x0, 1000, p); 

sol = solve(prob, KenCarp4(), reltol = 1e-12);