Getting Started

Installation

Run

add https://github.com/jo-ap/tikhonovfenichelreductions.jl

in Julia package Mode.

Example

Here we consider the derivation of the Rosenzweig-MacArthur model as a reduction from a three dimensional system as demonstrated in [8].

Load the package and its dependency Oscar.jl. Note that loading Oscar is optional, but results in prettier printing of types and imports useful functions.

using TikhonovFenichelReductions
using Oscar # optional

Define the components, parameters and the RHS of the system.

# components
x = ["B", "S", "H"]

# parameters
p = ["α", "β", "γ", "δ", "η", "ρ"]

# RHS of ODE system ẋ = f(x, p), where f is polynomial in x and p
function f(x, p)
  B, S, H = x
  α, β, γ, δ, η, ρ = p
  return [
    ρ*B*(1-B) - α*B*H,
    -η*S + γ*B*H,
    β*S- δ*H + η*S - γ*B*H
  ]
end

f (generic function with 1 method)

Initialize the problem with desired dimension of the reduced system

# dimension of the reduced system
s = 2

# create problem
problem = ReductionProblem(f, x, p, s)

ReductionProblem for dimension s = 2
 x = [B, S, H]
 p = [α, β, γ, δ, η, ρ]
 ODE system:
  dB/dt = -α*B*H - ρ*B^2 + ρ*B
  dS/dt = γ*B*H - η*S
  dH/dt = β*S - γ*B*H - δ*H + η*S

Slow-Fast Separations of Rates

We can find all slow-fast separations that are TFPVs by using necessary conditions for the existence of a reduction. This returns a vector of boolean indices for each TFPV candidate, where $0$ corresponds to a small parameters, i.e. those multiplied by $\varepsilon$, and $1$ corresponds to a parameter in $\mathcal{O}(1)$. In addition, we obtain the generators of the irreducible components of the affine variety

\[\mathcal{V}_{\mathbb{C}}(f(\cdot,\pi^\star)) = \{x\in\mathbb{C}^n \mid f(x,\pi^\star)=0\}\]

and their dimension, where $\pi^\star$ is defined by the corresponding slow-fast separation in sf_separation. This information is stored using the type Variety. Note that later have to check manually whether the variety taken in $\mathbb{R}^n$ has the same dimension (i.e. if there exists a real non-singular point), which then at least locally renders this variety a manifold.

tfpvs, varieties = tfpvs_and_varieties(problem)

(Vector{Bool}[[0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 1, 0], [0, 0, 0, 1, 0, 0], [0, 0, 1, 0, 0, 0], [0, 0, 1, 0, 0, 1], [0, 0, 1, 0, 1, 0], [0, 0, 1, 1, 0, 0], [0, 1, 0, 0, 1, 0], [0, 1, 0, 1, 0, 0], [1, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 1], [1, 0, 0, 1, 0, 0], [1, 0, 1, 0, 0, 0], [1, 0, 1, 0, 0, 1], [1, 0, 1, 1, 0, 0]], Vector{Variety}[[Variety(Ideal (B), Nemo.QQMPolyRingElem[B], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[B], Nemo.QQMPolyRingElem[B], [1], 2), Variety(Ideal (B - 1), Nemo.QQMPolyRingElem[B - 1], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[B - 1], Nemo.QQMPolyRingElem[B - 1], [1], 2)], [Variety(Ideal (S), Nemo.QQMPolyRingElem[S], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[S], Nemo.QQMPolyRingElem[S], [1], 2)], [Variety(Ideal (H), Nemo.QQMPolyRingElem[H], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[H], Nemo.QQMPolyRingElem[H], [1], 2)], [Variety(Ideal (H), Nemo.QQMPolyRingElem[H], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[H], Nemo.QQMPolyRingElem[H], [1], 2), Variety(Ideal (B), Nemo.QQMPolyRingElem[B], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[B], Nemo.QQMPolyRingElem[B], [1], 2)], [Variety(Ideal (B), Nemo.QQMPolyRingElem[B], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[B], Nemo.QQMPolyRingElem[B], [1], 2), Variety(Ideal (H, B - 1), Nemo.QQMPolyRingElem[H, B - 1], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[H, B - 1], Nemo.QQMPolyRingElem[H, B - 1], [1 0; 0 1], 1)], [Variety(Ideal (γ*B*H - η*S), Nemo.QQMPolyRingElem[γ*B*H - η*S], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[γ*B*H - η*S], Nemo.QQMPolyRingElem[γ*B*H - η*S], [1], 2)], [Variety(Ideal (H), Nemo.QQMPolyRingElem[H], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[H], Nemo.QQMPolyRingElem[H], [1], 2)], [Variety(Ideal (S), Nemo.QQMPolyRingElem[S], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[S], Nemo.QQMPolyRingElem[S], [1], 2)], [Variety(Ideal (β*S - δ*H), Nemo.QQMPolyRingElem[β*S - δ*H], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[β*S - δ*H], Nemo.QQMPolyRingElem[β*S - δ*H], [1], 2)], [Variety(Ideal (H), Nemo.QQMPolyRingElem[H], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[H], Nemo.QQMPolyRingElem[H], [1], 2), Variety(Ideal (B), Nemo.QQMPolyRingElem[B], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[B], Nemo.QQMPolyRingElem[B], [1], 2)], [Variety(Ideal (B), Nemo.QQMPolyRingElem[B], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[B], Nemo.QQMPolyRingElem[B], [1], 2), Variety(Ideal (ρ*B + α*H - ρ), Nemo.QQMPolyRingElem[α*H + ρ*B - ρ], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[ρ*B + α*H - ρ], Nemo.QQMPolyRingElem[α*H + ρ*B - ρ], [1], 2)], [Variety(Ideal (H), Nemo.QQMPolyRingElem[H], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[H], Nemo.QQMPolyRingElem[H], [1], 2)], [Variety(Ideal (H), Nemo.QQMPolyRingElem[H], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[H], Nemo.QQMPolyRingElem[H], [1], 2), Variety(Ideal (B), Nemo.QQMPolyRingElem[B], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[B], Nemo.QQMPolyRingElem[B], [1], 2)], [Variety(Ideal (B), Nemo.QQMPolyRingElem[B], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[B], Nemo.QQMPolyRingElem[B], [1], 2), Variety(Ideal (H, B - 1), Nemo.QQMPolyRingElem[H, B - 1], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[H, B - 1], Nemo.QQMPolyRingElem[H, B - 1], [1 0; 0 1], 1)], [Variety(Ideal (H), Nemo.QQMPolyRingElem[H], AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.RationalFunctionFieldElem{Nemo.QQFieldElem, Nemo.QQMPolyRingElem}}[H], Nemo.QQMPolyRingElem[H], [1], 2)]])

Show the results: All possible slow-fast separation of rates

print_tfpvs(problem, tfpvs)

π   = (α, β, γ, δ, η, ρ)
_______________________
π₁  = (⋅, ⋅, ⋅, ⋅, ⋅, ρ)
π₂  = (⋅, ⋅, ⋅, ⋅, η, ⋅)
π₃  = (⋅, ⋅, ⋅, δ, ⋅, ⋅)
π₄  = (⋅, ⋅, γ, ⋅, ⋅, ⋅)
π₅  = (⋅, ⋅, γ, ⋅, ⋅, ρ)
π₆  = (⋅, ⋅, γ, ⋅, η, ⋅)
π₇  = (⋅, ⋅, γ, δ, ⋅, ⋅)
π₈  = (⋅, β, ⋅, ⋅, η, ⋅)
π₉  = (⋅, β, ⋅, δ, ⋅, ⋅)
π₁₀ = (α, ⋅, ⋅, ⋅, ⋅, ⋅)
π₁₁ = (α, ⋅, ⋅, ⋅, ⋅, ρ)
π₁₂ = (α, ⋅, ⋅, δ, ⋅, ⋅)
π₁₃ = (α, ⋅, γ, ⋅, ⋅, ⋅)
π₁₄ = (α, ⋅, γ, ⋅, ⋅, ρ)
π₁₅ = (α, ⋅, γ, δ, ⋅, ⋅)

and the corresponding irreducible components containing the slow manifold and their Krull dimension

print_varieties(varieties)

V₁  : [B], 2
      [B - 1], 2
V₂  : [S], 2
V₃  : [H], 2
V₄  : [H], 2
      [B], 2
V₅  : [B], 2
      [H, B - 1], 1
V₆  : [γ*B*H - η*S], 2
V₇  : [H], 2
V₈  : [S], 2
V₉  : [β*S - δ*H], 2
V₁₀ : [H], 2
      [B], 2
V₁₁ : [B], 2
      [α*H + ρ*B - ρ], 2
V₁₂ : [H], 2
V₁₃ : [H], 2
      [B], 2
V₁₄ : [B], 2
      [H, B - 1], 1
V₁₅ : [H], 2

For further use, we need to make the variables (the system components and parameters) available in the Main namespace. Then, we can compute the reduction corresponding to TFPV 15, which is the Rosenzweig-MacArthur model.

B, S, H = system_components(problem)
α, β, γ, δ, η, ρ = system_parameters(problem)

# instantiate reduction
reduction = Reduction(problem, tfpvs[15])

Reduction for dimension s = 2 with
 slow: β, η, ρ
 fast: α, γ, δ

To find the slow manifold, we can consider the affine variety defined by the vanishing of

varieties[15][1].gens_R

1-element Vector{Nemo.QQMPolyRingElem}:
 H

which has (complex) dimension

varieties[15][1].dim

2

as a variety. We can see, that there is only one irreducible component, so the slow manifold is $\mathcal{V}_{\mathbb{R}}(H) = \{(B,S,0) \mid B,S \in \mathbb{R}\}$ with dimension $s=2$ as desired.

We need to define the slow manifold explicitly in order to check whether the reduction exists. This also allows us to substitute the variables that got reduced according to the slow manifold in the reduced system.

set_manifold!(reduction, [B, S, 0])

true

Note that in this case any generic point on the affine variety is non-singular and can be chosen. Thus, we don't have to call set_point! to set a non-singular point explicitly on whose neighbourhood the reduction exists. This also means $\dim \mathcal{V}_{\mathbb{C}}(f(\cdot,\pi^\star)) = \dim \mathcal{V}_{\mathbb{R}}(f(\cdot,\pi^\star))$ locally for the irreducible component containing the non-singular point.

Lastly, we define a product decomposition $f(\cdot,\pi^\star) = P \cdot \psi$ with $\mathcal{V}(\psi) = \mathcal{V}(f(\cdot,\pi^\star))$. This can be done by specifying only $\psi$ using the generators of the irreducible component of $\mathcal{V}(f^{(0)})$ that corresponds to the slow manifold. Then, the methods automatically computes $P$, which relies on having exactly $r$ generators of the variety corresponding to the slow manifold.

set_decomposition!(reduction, varieties[15][1])

true

In most cases, this method should work, but you can also manually set $P$ and $Psi$.

Now we can compute and show the reduced system.

compute_reduction!(reduction)
print_reduced_system(reduction)

dB/dt = -ρ*B^2 + ρ*B + (-α*β*B*S - α*η*B*S)//(γ*B + δ)
dS/dt = β*S + (-β*δ*S - δ*η*S)//(γ*B + δ)

This updates the reduction object, which now contains the reduced system before and after variables are substituted as defined by the slow manifold. You can see how the reduction object is updated:

reduction

Reduction for dimension s = 2 with
 slow: β, η, ρ
 fast: α, γ, δ
 M      = [B, S, 0]
 P      = [-α*B; γ*B; -γ*B-δ]
 Ψ      = [H]
 x₀     = [B, S, 0]
 Df(x₀) = [0 0 -α*B; 0 0 γ*B; 0 0 -γ*B-δ]
 Reduced system:
  dB/dt = -ρ*B^2 + ρ*B + (-α*β*B*S - α*η*B*S)//(γ*B + δ)
  dS/dt = β*S + (-β*δ*S - δ*η*S)//(γ*B + δ)

The first two components of reduction.g define the reduced system and can be rewritten as

\[\begin{align*} \frac{dB}{dt} &= \rho B (1 - B) - \alpha(\eta + \beta) S \frac{B}{\delta + \gamma B} \\ \frac{dS}{dt} &= -\eta S + \gamma(\eta + \beta) S \frac{B}{\delta + \gamma B} \\ \end{align*}\]

which is exactly the Rosenzweig-MacArthur model.

The slow manifold is attractive if all non-zero eigenvalues of the Jacobian of $f(\cdot, \pi^\star)$ at the non-singular x0 point have negative real part.

# pretty-print Jacobian at x0
J = jacobian_tfpv_at_x0(reduction)
show(stdout, "text/plain", J)

[0   0       -α*B]
[0   0        γ*B]
[0   0   -γ*B - δ]

Thus, the full system converges to the reduction as $\varepsilon \to 0$ if $B \geq 0$.

Bulk Computations

TikhonovFenichelReductions.jl has methods that simplify the computation of multiple reductions with the same manifold, since typically many different TFPVs share the same slow manifold. Additionally, the function get_explicit_manifold implements a heuristic to find an explicit parametric description of the slow manifold. Thus, in simple enough cases, the computation of all reductions is fully automatic.

all_V = unique_varieties(problem, varieties)
all_M = [get_explicit_manifold(problem, V) for V in all_V]

# get reduction and indices of varieties that correspond to unique slow manifolds
all_reductions, idx_M = compute_all_reductions(problem, tfpvs, varieties, [m[1] for m in all_M]);

all_reductions[15][1]

Reduction for dimension s = 2 with
 slow: β, η, ρ
 fast: α, γ, δ
 M      = [B, S, 0]
 P      = [-α*B; γ*B; -γ*B-δ]
 Ψ      = [H]
 x₀     = [B, S, 0]
 Df(x₀) = [0 0 -α*B; 0 0 γ*B; 0 0 -γ*B-δ]
 Reduced system:
  dB/dt = -ρ*B^2 + ρ*B + (-α*β*B*S - α*η*B*S)//(γ*B + δ)
  dS/dt = β*S + (-β*δ*S - δ*η*S)//(γ*B + δ)

General TFPVs

One can also get all general TFPVs by computing a Gröbner basis G that reflects necessary conditions on the parameters of the system. Note that depending on the input system and the drop in dimension this can be a very computationally intensive task.

G = tfpvs_groebner(problem)

4-element Vector{Nemo.QQMPolyRingElem}:
 δ*η*ρ
 γ*δ*ρ
 β*γ*ρ
 α*δ*η

The TFPVs are points in the affine variety of G. Thus, in cases where G is simple enough, we can compute a minimal primary decomposition to characterise the cases in which G vanishes. Note that here we use functions exported by Oscar.jl.

I = ideal(G)
PD = primary_decomposition(I)
for Q in PD
    println(gens(Q[2]))
end

Nemo.QQMPolyRingElem[ρ, η]
Nemo.QQMPolyRingElem[ρ, δ]
Nemo.QQMPolyRingElem[ρ, α]
Nemo.QQMPolyRingElem[η, γ]
Nemo.QQMPolyRingElem[δ, γ]
Nemo.QQMPolyRingElem[δ, β]

Thus, in this case, all TFPVs are slow-fast separations of rates.

In cases where G is more complicated, we can check whether the ideal generated by G contains any monomials as follows.

# polynomial ring ℚ[p,x]
R = parent(α)
# polynomial ring ℚ[p]
S, v = polynomial_ring(QQ, "_" .* p)
h = hom(S, R, system_parameters(problem))

# the ideal generated by G in the ring S
I = preimage(h, ideal(G));

# compute the saturation I:⟨π₁⋯πₘ⟩^∞
I_sat = saturation(I, ideal(prod(v)))
is_one(I_sat)

true