Change parameters in a model in specific time with DifferentialEquations

[askery]

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:

DifferentialEquations with nested functions + Callbacks

Hello, now that I (sort of) know how to calculate differential equations with nested functions I would like to learn how to use it with callbacks. Callbacks are used to run a differential equation to a certain point, changing the value of a parameter and then continue running. And that’s really handy when working with seasonal data. I know how a callback looks like when using a differential equation that does not have nested functions, but my experiments with a callback and differential equatio…

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

Time-dependent events in ODE

You need to split into two callbacks (as your code errors on this), see http://docs.juliadiffeq.org/latest/features/callback_functions.html#CallbackSet-1 on how to combine them. Then the other error is function affect!(integrator) integrator.u[1] = stim end as stated in the docs: http://docs.juliadiffeq.org/latest/features/callback_functions.html#Example-1:-Bouncing-Ball-1 So all in all, slightly modified: using Plots using DifferentialEquations # Parameters const k21 = 0.14*24 const k…

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.

[ChrisRackauckas]

@askery

a = 2

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

[askery]

@askery

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)

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)

[eliassno]

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

---

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)

[askery]

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

[peterj]

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

[askery]

Cool, will try that too. Thanks, peterj.

[fabern]

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]

[fabern]

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?