Generating zero-isoclines for a coupled ODE

General Usage
1

Hi All,

This seems like it should be a simple problem to solve, but I’ve scoured the documentation of DifferentialEquations.jl and can’t seem to identify anything specifically relating to my issue.

For a simple coupled ODE of the form

\frac{dA}{dt} = f(A,B,t)\\ \frac{dB}{dt} = g(A,B,t)

I’ve set up my Julia implementation like so:

u0 = [20.0, 2.6]; #Initial populations
p = [ib, r, nutr, H₁, H₂, H₃, H₄];
tspan = (0.0,364.0);
prob = ODEProblem(coupled_ode!,u0,tspan,p);
sol = solve(prob);

I want to generate a state space diagram over A \in [0,20], B \in [0,120], showing zero isoclines for both parameters (i.e. traces over the state space when \frac{dA}{dt} = 0 or \frac{dB}{dt} = 0).

There are a few steady state solutions to this problem, that I can obtain using prob = SteadyStateProblem(coupled_ode!,u0,p) in place of the ODEProblem, but I can’t seem to find a concise/efficient solution for generating the zero-isoclines at all.

Currently I’ve looked into using two for loops altering u0 = [A0[i], B0[j]] to build the state space.
I can use sol(0.0, Val{1}) for each of the solutions to find the derivative values when t=0, but of course I’m looking for the the time when Val{1} == 0.
Perhaps I can generate a large search range over time and find the index of the minimum of the derivative?

search = 0.0:0.001:5000.0;
findmin(abs.(hcat(sol(search,Val{1}).u...)),dims=2)

but this seems like overkill, and doesn’t necessarily assure that A’ = 0 has been found either.

It’s highly likely I’ve either missed a simple method or misunderstood something conceptually. Any hints or suggestions on this would be appreciated.

2

@rveltz I assume there’s a nice way to do this with [ PseudoArcLengthContinuation.jl](https://rveltz.github.io/PseudoArcLengthContinuation.jl/dev/) ?

2 Likes
3

you can try

using PseudoArcLengthContinuation
const PALC = PseudoArcLengthContinuation

h = (A,B) -> [A^2-B]
contparams = ContinuationPar(pMin = -1., pMax = 1.,
                ds = 0.001) # set ds = -0.001 for the other part of the manifold
nullcline, _,_ = continuation( (x,p)->h(x[1],p),
                [0.], #value of A for B = 0
                0., # initial value of B
                contparams; printSolution = (x, p) -> x[1])

plot(nullcline.branch[2,:], nullcline.branch[1,:])

Another idea: https://github.com/bachrathyd/MDBM.jl

1 Like
4

If f and g are independent of t (which they must be, otherwise I don’t think the concept of zero-isocline makes sense?) then you are just asking to draw contours or level sets of the functions f and g:

using Plots
f(x, y) = x^2 + 2y^2
g(x, y) = x - 2y^3

myrange = -2:0.01:2

contour(myrange, myrange, f, levels=[0])
contour!(myrange, myrange, g, levels=[0])
5 Likes
5

Thanks everyone.

@dpsanders, you are correct. It’s my conceptual issue with this problem that was the killer. Thanks!

@ChrisRackauckas & @rveltz: will be interesting to take a look at the packages you suggested for other problems I have in the future - thank you. This one however didn’t even need diffeq in the end (other than for identification of the solutions to the SteadyStateProblem().

1 Like