Help finding multiple roots to nonlinear equation - IntervalRootFinding not working

[lpz18]

I have a set of nonlinear equations (which happen to be steady-state DEs) that I am trying to solve that have multiple real roots and 3 variables (x, y, & z in the code below)

In Mathematica I can get the multiple roots as I’d expect:

In Julia, using JuMP I cannot:

(the difference I care about in comparing the two graphs is that the curves go jagged though they shouldn’t. the axes are slightly different since I re-scaled in the Julia figure, and the points are data which can be neglected.)

I would like to find the multiple solutions I have in Mathematica. Naturally, JuMP finds the single optimal according to the specified objective, so this is as expected but not as desired.

I tried IntervalRootFinding.jl:

function core_eqns(u)
	x, y, z = u

	φ_R = 2 * 0.30 * x / (x + 0.8) + 0.05
	φ_F = 0.35 - φ_R
	φ_M = φ_F * 0.99
	φ_AR = φ_F * 0.01
	C_Rf = 63.7e-6 * φ_R - y
	φ_Rf = φ_R - y / 63.7e-6
	λ = 10 * x / (x + 0.03) * φ_Rf

	# steady-state equations
	A = 5 * φ_M - (1 + x) * λ
	B = 2.25e6 * z * C_Rf - 3.78 * y - λ * y
	C = 800 * (5e-4 - z) - 35 * 60^2 * φ_AR * 536e-6 * z / (z + 15e-6) + 3.78 * y - 2.25e6 * C_Rf * z - λ * z
	
	return [A, B, C]
end

function SS_solver(p)
	A_range = 0..15
	B_range = 0..35e-6
	C_range = 0..5e-4
	interval = A_range × B_range × C_range

	function f(x)
		out = core_eqns(x, p)
		return SVector(
			out[1],
			out[2],
			out[3]
		)
	end
	return roots(f, interval)
end

this gave me ~2000 intervals all with the flag :unknown, and it took 455s to run – which is quite undesirable since each the curves above are generated from computing the roots 100 times (and for comparison it takes <4s for all of the solutions to run in JuMP)

[kapple]

What’s the result from Roots.jl?

[lpz18]

from the docs:

This package contains simple routines for finding roots, or zeros, of scalar functions of a single real variable

I didn’t try because I have 3 variables. Is there a way of using Roots.jl for multivariate root-finding?

[kapple]

Ah I guess not then.
Sorry, I’m on mobile and just messaged a quick suggestion.

[lpz18]

thanks anyways

[kapple]

How about this?

Steady State Callbacks · DiffEqCallbacks.jl

[lpz18]

the DE steady solvers are the best combination of speed and robustness that I’ve found, but haven’t been able to get the multiple real solutions, I might just need to more systematically sweep different initial conditions

[orebas]

If these are rational, you might try clearing denominators and then HomotopyContinuation.jl.

[kapple]

@ChrisRackauckas any ideas?

[ChrisRackauckas]

Isn’t this just a bifurcation plot? Why not use BifurcationKit.jl?

[lpz18]

HomotopyContinuation.jl and BifurcationKit.jl both look promising! I’ll check them out and report back, thanks for the help!

[lpz18]

With BifurcationKit.jl I have not managed to find a bifurcation, this is my code:

	φ_F = 1 - 0.65 - φ_R
	φ_M = φ_F * 0.99
	φ_AR = φ_F * 0.01
	C_Rf = 63.7e-6 * φ_R - C_bnd
	φ_Rf = φ_R - C_bnd / 63.7e-6
	λ = 9.65 * c_pc / (c_pc + 0.03) * φ_Rf

	# steady-state equations
	cpc_ss = 5 * φ_M - (1 + c_pc) * λ
	C_bnd_ss = 2.25e6 * C_in * C_Rf - 3.78 * C_bnd - λ * C_bnd
	C_in_ss = 760 * (abx_conc[1] - C_in) - 35*60^2 * φ_AR * 536e6 * C_in / (C_in + 15e-6) + 3.78 * C_bnd - 2.25e6 * C_Rf * C_in - λ * C_in
	
	return [cpc_ss, C_bnd_ss, C_in_ss]
end

prob = BifurcationProblem(core_eqns, zeros(3), [1e-5], 1)
br = continuation(prob, PALC(), ContinuationPar(
	p_min = 0.0, p_max = 1e-3, detect_bifurcation = 2, max_steps = 1000
))

and the output I get is:

 ┌─ Curve type: EquilibriumCont
 ├─ Number of points: 2
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter p starts at 1.0e-5, ends at 0.001
 ├─ Algo: PALC
 └─ Special points:

- #  1, endpoint at p ≈ +0.00100000,   

[lpz18]

HomotopyContinuation.jl is working for single solutions but I am still working out how to run the loop so I can plot the multiple solutions nicely like in the original question. Strictly speaking, HomotopyContinuation.jl has satisfied my question as originally stated, but I still need to get the analogous plot which requires looping through the solution multiple times (like is intended in the bifurcation analysis)

[lpz18]

got it working nicely with BifurcationKit.jl, thanks for the pointers