# Switching behaviour of implicit ODEs

**URL:** https://discourse.julialang.org/t/switching-behaviour-of-implicit-odes/121472
**Category:** Modelling & Simulations
**Created:** [October 18, 2024, 11:28pm UTC](https://discourse.julialang.org/t/switching-behaviour-of-implicit-odes/121472 "2024-10-18T23:28:20Z")
**Posts on this page:** 8
**Page:** 1

<div class="post-metadata">

### Author: ![clm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/clm/32/45923_2.png) [@clm](https://discourse.julialang.org/u/clm)
#### Post date: [October 18, 2024, 11:28pm UTC](https://discourse.julialang.org/t/switching-behaviour-of-implicit-odes/121472/1 "2024-10-18T23:28:20Z")

</div>

I have two coupled nonlinear ODEs where I’m aiming to find x(t) and y(t), defined generally as:  
f(x,y,x',y') = 0,  
\min(g(x,y,x',y'),h(x,y,x',y')) = 0,

where I expect that I will switch between g\<h and g\>h at least once at some unknown location in the domain.

Naïvely, I first tried the following (forgive the clunky variable names)

```julia
using DifferentialEquations, Plots, Sundials
function test_implicit(res,du,u,p,z)
    # unpack for visual clarity
    κ, Ḣ = p
    θ̃,τ̃ = u
    θ̃_z,τ̃_z = du
    
    func_g = Ḣ*θ̃ + κ*θ̃_z ### g(...)
    func_h = Ḣ^2*θ̃^2 - τ̃*θ̃_z### h(...)
    
    res[1] = τ̃*θ̃_z - τ̃_z*θ̃ - θ̃^2 ### f(...) = 0
    res[2] = min(condA,condB)
    
end

κh = 2.5#1.8
Ḣ0 = 5.0
p = κ, Ḣ

zmin = 0.
zmax = 3.

prob = DAEProblem(test_implicit,[0.0,0.0], [-1.0,4.0], (zmin, zmax),
        p, differential_vars = [true,true], saveat=0.01)

sol = solve(prob1,DFBDF());
plot(reduce(hcat,sol.u)[1,:],sol.t)
plot!(reduce(hcat,sol.u)[2,:],sol.t)

```

At first, it appears to work – but that’s because the initial parameters \kappa and \dot{H} happen to define a regime where g\<h everywhere. If instead we set \kappa=0.5, the integrator doesn’t return anything, as a ``regime change" occurs in the domain. We can confirm this by setting `res[2] = func_g` and `res[2] = func_h` in separate runs, and see that both run, and neither case blows up.

At this point I’m pretty certain I need to be using ContinuousCallbacks, but all of the examples I found were for explicit, not implicit function definitions. So taking the bouncing ball example code (just to make something work, it won’t do anything meaningful here),

```julia
condition(out, du, u, p, z, integrator) = u[1]
function affect!(integrator)
    integrator.u[2] = -integrator.u[2]
end

cb = ContinuousCallback(condition, affect!)

```

I get the error,

```julia
MethodError: no method matching test_implicit(::Vector{Float64}, ::Vector{Float64}, ::Tuple{Float64, Float64}, ::Float64)
Closest candidates are:
  test_implicit(::Any, ::Any, ::Any, ::Any, ::Any)

```

It seems I can’t supply enough arguments to my condition function. I admit I don’t really understand the integrator interface, and I also notice that none of the examples with ContinuousCallback use du in the `condition` (which I will need to do, as g and h depend on the derivatives), so I hope that is possible.

Thanks in advance!

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 19, 2024, 12:47am UTC](https://discourse.julialang.org/t/switching-behaviour-of-implicit-odes/121472/2 "2024-10-19T00:47:43Z")

</div>

> [@clm](#):
>
> `condition(out, du, u, p, z, integrator) = u[1]`

You added extra arguments that don’t exist in there. `condition(u, z, integrator) = u[1]`

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 19, 2024, 12:48am UTC](https://discourse.julialang.org/t/switching-behaviour-of-implicit-odes/121472/3 "2024-10-19T00:48:03Z")

</div>

> [@clm](#):
>
> ContinuousCallback use dududu in the `condition` (which I will need to do, as ggg and hhh depend on the derivatives), so I hope that is possible.

`integrator(t, Val{1})`

---

<div class="post-metadata">

### Author: ![clm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/clm/32/45923_2.png) [@clm](https://discourse.julialang.org/u/clm)
#### Post date: [October 19, 2024, 5:48am UTC](https://discourse.julialang.org/t/switching-behaviour-of-implicit-odes/121472/4 "2024-10-19T05:48:40Z")

</div>

Sorry, I should have made this more clear. I had tried this and it wasn’t working initially (hence the attempt at incorrect arguments).  
I now find that, say,  
`function condition(u, z, integrator) = u[1]-0.5`  
works fine (this condition condition that is never met, as u[1] goes from -1 to 0). But,

`function condition(u, z, integrator) = u[1]+0.5`

Throws the same error as above, when using an example affect,

```julia
function affect!(integrator)
    integrator.u[1] = -integrator.u[1]
end

```

So maybe my issue is in the affect function?

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 19, 2024, 6:06am UTC](https://discourse.julialang.org/t/switching-behaviour-of-implicit-odes/121472/5 "2024-10-19T06:06:52Z")

</div>

> [@clm](#):
>
> ```julia
> MethodError: no method matching test_implicit(::Vector{Float64}, ::Vector{Float64}, ::Tuple{Float64, Float64}, ::Float64)
> 
> ```

That error is completely unrelated.

Did you try it with IDA?

---

<div class="post-metadata">

### Author: ![clm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/clm/32/45923_2.png) [@clm](https://discourse.julialang.org/u/clm)
#### Post date: [October 19, 2024, 7:45pm UTC](https://discourse.julialang.org/t/switching-behaviour-of-implicit-odes/121472/6 "2024-10-19T19:45:29Z")

</div>

Thanks, Chris. I thought the error seemed unrelated but didn’t identify that the integrator choice was the cause. IDA and `integrator(t, Val{1})` was the ticket. I had even read that section of the integrator interface docs, but didn’t quite piece it together.

I’ll put my final code here for posterity in case anyone else comes across this thread. I’m using a ContinuousVectorCallback to check for the conditions on differentiable\_min(A,B)-A and differentiable\_min(A,B)-B and using that to update a switching parameter between 0 and 1. Requires a small reltol and abstol to run nicely.

```julia
function differentiable_min(A,B)
    (A+B)/2 - sqrt((A-B)^2)/2
end
    
function test_implicit(res,du,u,p,z)
    κh, Ḣ0, p_switch = p
    θ̃,τ̃ = u
    θ̃_z,τ̃_z = du
            
    condA = Ḣ0*θ̃ + κh*θ̃_z
    condB = Ḣ0^2*θ̃^2 - τ̃*θ̃_z
    
    
    res[1] = condA * p_switch + condB * (1 - p_switch)
    res[2] = τ̃*θ̃_z - τ̃_z*θ̃ - θ̃^2
end

function condition(out, u, z, integrator)    
    condA = integrator.p[2]*integrator.u[1] + integrator.p[1]*integrator(z, Val{1})[1]
    condB = integrator.p[2]^2*integrator.u[1]^2 - integrator.u[2]*integrator(z, Val{1})[1]
    
    out[1] = differentiable_min(condA,condB)-condA
    out[2] = differentiable_min(condA,condB)-condB
end

function affect!(integrator, idx)
    if idx == 1
        integrator.p[3] = 1
    elseif idx == 2
        integrator.p[3] = 0
    end
end

cb = VectorContinuousCallback(condition, affect!, 2)

κh = 0.3#1.8
Ḣ0 = 5.0
regime_init = 1
p = [κh, Ḣ0, regime_init]

zmin = 0.
zmax = 0.5

prob1 = DAEProblem(test_implicit,[1.0,-1.0], [-1.0,4.0], (zmin, zmax),
        p, differential_vars = [true,true])

```

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 19, 2024, 9:38pm UTC](https://discourse.julialang.org/t/switching-behaviour-of-implicit-odes/121472/7 "2024-10-19T21:38:23Z")

</div>

> [@clm](#):
>
> I thought the error seemed unrelated but didn’t identify that the integrator choice was the cause.

DFBDF still has an issue open about some parts of its event handling. I should get to that soon, but since the fastest thing is always to use the mass matrix form it hasn’t been a priority.

---

<div class="post-metadata">

### Author: ![clm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/clm/32/45923_2.png) [@clm](https://discourse.julialang.org/u/clm)
#### Post date: [October 20, 2024, 12:19am UTC](https://discourse.julialang.org/t/switching-behaviour-of-implicit-odes/121472/8 "2024-10-20T00:19:51Z")

</div>

I thought I had this sorted, but it seems my solution is highly susceptible to my choice of dt. Changing the reltol and abstol from `1e-13` to `1e-12` changes all of these curves slightly, which makes me think none of them are the “true” solution.

![tau_curves](https://global.discourse-cdn.com/julialang/original/3X/c/0/c076d78da7e04bf4e6519fc1e4f0fe181c953ecd.png)

Figure inelegantly made with:

```julia
p = [κh, Ḣ0, regime_init]
prob1 = DAEProblem(test_implicit,[1.0,-1.0], [-1.0,4.0], (zmin, zmax),
        p, differential_vars = [true,true])
sol = solve(prob1,IDA(),callback=cb,abstol=1e-13,reltol=1e-13);
plot(reduce(hcat,sol.u)[2,:],sol.t,label="dz determined automatically",c="pink",ylabel="z",xlabel=L"$\tilde{\tau}$",linewidth=5,xlim=(0,1))

p = [κh, Ḣ0, regime_init]
prob1 = DAEProblem(test_implicit,[1.0,-1.0], [-1.0,4.0], (zmin, zmax),
        p, differential_vars = [true,true],dt=0.0001)
sol = solve(prob1,IDA(),callback=cb,abstol=1e-13,reltol=1e-13);
plot!(reduce(hcat,sol.u)[2,:],sol.t,label="dz=0.0001",ylabel="z",c="blue",xlabel=L"$\tilde{\tau}$",linewidth=5,xlim=(0,1))

p = [κh, Ḣ0, regime_init]
prob1 = DAEProblem(test_implicit,[1.0,-1.0], [-1.0,4.0], (zmin, zmax),
        p, differential_vars = [true,true],dt=0.00001)
sol = solve(prob1,IDA(),callback=cb,abstol=1e-13,reltol=1e-13);
plot!(reduce(hcat,sol.u)[2,:],sol.t,label="dz=0.00001",linewidth=3)

p = [κh, Ḣ0, regime_init]
prob1 = DAEProblem(test_implicit,[1.0,-1.0], [-1.0,4.0], (zmin, zmax),
        p, differential_vars = [true,true],dt=0.000001)
sol = solve(prob1,IDA(),callback=cb,abstol=1e-13,reltol=1e-13);
plot!(reduce(hcat,sol.u)[2,:],sol.t,label="dz=0.000001",linewidth=3, c="orange")

p = [κh, Ḣ0, regime_init]
prob1 = DAEProblem(test_implicit,[1.0,-1.0], [-1.0,4.0], (zmin, zmax),
        p, differential_vars = [true,true],dt=0.0000001)
sol = solve(prob1,IDA(),callback=cb,abstol=1e-13,reltol=1e-13);
plot!(reduce(hcat,sol.u)[2,:],sol.t,label="dz=0.0000001",linewidth=3,ls=:dash,c=:black)

p = [κh, Ḣ0, regime_init]
prob1 = DAEProblem(test_implicit,[1.0,-1.0], [-1.0,4.0], (zmin, zmax),
        p, differential_vars = [true,true],dt=0.00000001)
sol = solve(prob1,IDA(),callback=cb,abstol=1e-13,reltol=1e-13);
plot!(reduce(hcat,sol.u)[2,:],sol.t,label="dz=0.00000001",linewidth=3,ls=:dash,c=:grey)

p = [κh, Ḣ0, regime_init]
prob1 = DAEProblem(test_implicit,[1.0,-1.0], [-1.0,4.0], (zmin, zmax),
        p, differential_vars = [true,true],dt=0.000000001)
sol = solve(prob1,IDA(),callback=cb,abstol=1e-13,reltol=1e-13);
plot!(reduce(hcat,sol.u)[2,:],sol.t,label="dz=0.000000001",linewidth=3,ls=:dash,c=:purple)

p = [κh, Ḣ0, regime_init]
prob1 = DAEProblem(test_implicit,[1.0,-1.0], [-1.0,4.0], (zmin, zmax),
        p, differential_vars = [true,true],dt=0.00000000001)
sol = solve(prob1,IDA(),callback=cb,abstol=1e-13,reltol=1e-13);
plot!(reduce(hcat,sol.u)[2,:],sol.t,label="dz=0.00000000001",linewidth=2,ls=:dashdot,c=:red)

```