# 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 = a*(u-u)
du = u*(b-u) - u
du = u*u - c*u
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 = a*(u-u)
du = u*(b-u) - u
du = u*u - c*u
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?

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

1 Like

Thanks, Chris for the prompt reply. Tried:

``````function lorenz!(du,u,p,t)
du = a*(u-u)
du = u*(b-u) - u
du = u*u - c*u
end
``````

And

``````function affect!(integrator)
integrator.p = 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 = a*(u-u)
du = u*(b-u) - u
du = u*u - c*u
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 = a*(u-u)
du = u*(b-u) - u
du = u*u - c*u
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
println("t = \$t, \t a = \$a")
end
return t < 50
end

function affect!(integrator)
integrator.p = 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)
``````
5 Likes

Thank you very much, eliassno. This is what I need 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.

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 = a*(u-u)
du = u*(b-u) - u
du = u*u - c*u
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 = 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?