Switching behaviour of implicit ODEs

Specific Domains Modelling & Simulations
1

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)

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),

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,

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!

2

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

3

integrator(t, Val{1})

4

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,

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

So maybe my issue is in the affect function?

5

That error is completely unrelated.

Did you try it with IDA?

6

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.

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])
7

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.

8

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

Figure inelegantly made with:

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)