Change parameters in a model in specific time with DifferentialEquations

Guys, I am really new to Julia. I came to solve a system of ODE and want to stay. For sure my question is trivial, but I looked for similar solutions and could not find for 3 days.

Suppose I am solving Lorentz equations. But at a given time I need to change a parameter:
say a = 10 to a = 3 at time t = 50.

How can I do that? I know I can use callbacks, but could not find minimal working examples.

using Plots
using DifferentialEquations

# parameters
a = 10.0
b = 28.0
c = 8/3

function lorenz!(du,u,t)
 du[1] = a*(u[2]-u[1])
 du[2] = u[1]*(b-u[3]) - u[2]
 du[3] = u[1]*u[2] - c*u[3]
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)
@time sol = solve(prob)

plot(sol,vars=(0,1))

Once again, I need to change a=10 to a=2, for example, in a given time t = 50.

Tried:

This was close to solve, where one can change the variables but not the parameters.

Inspired by the solution above I tried:

using Plots
using DifferentialEquations

# parameters
a = 10.0
b = 28.0
c = 8/3

function lorenz!(du,u,t)
 du[1] = a*(u[2]-u[1])
 du[2] = u[1]*(b-u[3]) - u[2]
 du[3] = u[1]*u[2] - c*u[3]
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)

# Events
function condition(y, t, integrator)
    t - 50
end
function affect!(integrator)
    a = 2
end
cb = ContinuousCallback(condition, affect!)
#cb = DiscreteCallback(condition, affect!)


# Solve
@time sol = solve(prob, Rodas4(), callback = cb)

I expected to see a = 2, but it does not change, it keeps a = 10.

Could someone help me?

Thanks in advance.

You need to have p in your differential equation, and then do integrator.p[1] = 2

1 Like

Thanks, Chris for the prompt reply. Tried:

function lorenz!(du,u,p,t)
 du[1] = a*(u[2]-u[1])
 du[2] = u[1]*(b-u[3]) - u[2]
 du[3] = u[1]*u[2] - c*u[3]
end

And

function affect!(integrator)
    integrator.p[1] = 2
end

But got this error (I am using julia 1.0.4 btw):

MethodError: no method matching setindex!(::DiffEqBase.NullParameters, ::Int32, ::Int32)

In lorenz!(du,u,p,t), you’ll want to unpack the parameters.

function lorenz!(du,u,p,t)
 a,b,c = p
 du[1] = a*(u[2]-u[1])
 du[2] = u[1]*(b-u[3]) - u[2]
 du[3] = u[1]*u[2] - c*u[3]
end

Then, you need to specify the parameters and add them to the ODEProblem.

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
p = [a,b,c]
prob = ODEProblem(lorenz!,u0,tspan,p)

I assume you want to set a = 10.0 when t < 50 and a = 2.0 when t > 50,
In this case, you should also change condition(y, t, integrator):

function condition(y, t, integrator)
     t < 50
end

I get the following
without_callback

with_callback


Below you can see the full code and output.

Full code
using Plots
using DifferentialEquations

# parameters
a = 10.0
b = 28.0
c = 8/3

function lorenz!(du,u,p,t)
 a,b,c = p
 du[1] = a*(u[2]-u[1])
 du[2] = u[1]*(b-u[3]) - u[2]
 du[3] = u[1]*u[2] - c*u[3]
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
p = [a,b,c]
prob = ODEProblem(lorenz!,u0,tspan,p)

# Events
function condition(y, t, integrator)
    if t < 50.03 && t > 49.96
        a = integrator.p[1]
        println("t = $t, \t a = $a")
    end
    return t < 50
end

function affect!(integrator)
    integrator.p[1] = 2.0
end
cb = ContinuousCallback(condition, affect!)


# Solve
@time sol = solve(prob, Rodas4())
plot(sol, title = "Without callback")
savefig("without_callback")


@time sol = solve(prob, Rodas4(), callback = cb)
plot(sol, title = "With callback")
savefig("with_callback")
Output
  2.182867 seconds (3.44 M allocations: 144.341 MiB, 2.40% gc time)
t = 49.969386198264516,          a = 10.0
t = 49.9618451821885,    a = 10.0
t = 49.96561569022651,   a = 10.0
t = 49.969386198264516,          a = 10.0
t = 49.969386198264516,          a = 10.0
t = 50.00655075541622,   a = 10.0
t = 50.00655075541622,   a = 10.0
t = 50.00655075541622,   a = 2.0
t = 50.00655497964498,   a = 2.0
t = 50.01124434293082,   a = 2.0
t = 50.01593793044543,   a = 2.0
t = 50.020631517960034,          a = 2.0
t = 50.025325105474636,          a = 2.0
  1.602665 seconds (2.59 M allocations: 112.577 MiB, 1.55% gc time)
3 Likes

Thank you very much, eliassno. This is what I need :grinning:

You could also use a conditional statement in lorenz! that checks the time and assigns a value for a accordingly.

2 Likes

Cool, will try that too. Thanks, peterj.

Instead of using ContinuousCallback, which may be called multiple times, you can get your desired result in a more straightforward way by using PresetTimeCallback which should be called only (and precisely) at t=50.0.

You can read about different callback types in: https://github.com/SciML/DiffEqCallbacks.jl#presettimecallback

Please see the MWE below.

MWE
using Plots
using DifferentialEquations

# parameters
a = 10.0
b = 28.0
c = 8/3

function lorenz!(du,u,p,t)
 a,b,c = p
 du[1] = a*(u[2]-u[1])
 du[2] = u[1]*(b-u[3]) - u[2]
 du[3] = u[1]*u[2] - c*u[3]
end

u0 = [1.0;0.0;0.0];
tspan = (0.0,100.0);
p = [a,b,c];
prob = ODEProblem(lorenz!,u0,tspan,p)

# Events variant 2
function affect!(integrator)
    integrator.p[1] = 2.0;
end
cb_variant2 = PresetTimeCallback(50.,affect!);

# Solve
print("\np before solving: $p\n")
print("\nprob.p before solving: $p\n")

@time sol = solve(prob, Rodas4(), callback = cb_variant2)
plot(sol, title = "With PresetTimeCallback")

print("\np after solving: $p\n")
print("\nprob.p after solving: $p\n")
MWE output
p before solving: [10.0, 28.0, 2.6666666666666665]

prob.p before solving: [10.0, 28.0, 2.6666666666666665]
2.090013 seconds (3.74 M allocations: 159.063 MiB, 2.91% gc time)

p after solving: [2.0, 28.0, 2.6666666666666665]

prob.p after solving: [2.0, 28.0, 2.6666666666666665]

Regarding above MWE I have myself a follow up question. affect! seems to modify the objects p and prob.p in the global scope. A subsequent run with those would thus use the modified value of a=2.0 from the start instead of the initial value.

Is there a way to modify a parameter value using a callback only for the current solver run so that a rerun is possible?