Differential-algebraic equations (DAEs) play an important role in dynamical system modeling such as power systems (cf. [1, 2]), nonlinear-circuits [3–7], robotics (cf. [8]), flight control systems [9], multi-body systems [10], numeric PDEs ([11] and references in Kunkel and Mehrmann [12]).

We consider quasilinear DAEs of form (1)A(x,α)x.=f(x,α), (x,α)∈ℝn×ℝm, for smooth functions A : ℝⁿ × ℝ^m → ℝ^{n × n}, f : ℝⁿ × ℝ^m → ℝⁿ of the phase variable x ∈ ℝⁿ and the parameter α ∈ ℝ^m.

In the presence of singularities, that is, points (x, α) such that det A(x, α) in Equation (1) vanishes, it is not possible to describe the local behavior of a DAE in terms of an explicit ODE. Regularization of such a singular DAE often leads to an ODE with higher dimensional manifolds of equilibria in phase space, which can manifest bifurcations without parameters (cf. [13]). System with singularities possessing the Hamiltonian structure were studied in [14–18].

In parametrized problems, a stability change due to the divergence of an eigenvalue was first analyzed by Venkatasubramanian [19, 20] and later addressed by many others [9, 20–25]. The change of stability, termed singularity-induced bifurcations (SIB), occurs when an equilibrium branch intersects a singular manifold, which results in the divergence of at least one eigenvalue through infinity.

Main efforts in studying such singularity-crossing equilibria have been given by trying to characterize the SIBs in terms of the linearized problems, such as using the matrix pencils {A(x*,α*),-Dxf(x*,α*)} associated to Equation (1) at the point of singularity (x^*, α^*) (cf. [20, 26]). Different sufficient conditions have been given in the framework of the tractability index (cf. [27, 28]) and the geometric index [26, 29, 30] among others. However, they may not provide a necessary and sufficient characterization of the local flow around such singularity-crossing sequilibria.

Besides SIBs, there can be other singular behavior induced by the presence of singularities such as the change of the singularity surface itself (fold) or bifurcations of singular equilibria that change significantly the dynamics near the singularity.

Quasilinear DAEs (Equation 1) have a strong connection to another important class of dynamical systems–fast-slow systems. Indeed, a system εx.=f(x,y), y.=g(x,y) is a fast-slow system for small ε. Setting ε to zero we obtain the so-called slow system: (2)0=f(x,y), y.=g(x,y), in which the first equation defines the slow manifold and the second one determines the dynamics, restricted onto it. System (2) is a DAE that can be brought to the form of an ODE system via time-differentiation of the algebraic equation and the substitution of ẏ from the second one: (3)fx(x,y)x.=-fy(x,y)g(x,y), y.=g(x,y). That is, the slow system can be reduced to the form (Equation 1), so our research here contributes also to the studies of bifurcations in slow-fast systems (see [31]). The correspondence between the terms and notions can be viewed in the following way:

A singularity set in DAEs corresponds to a fold set of a slow manifold;

In this paper, we will call bifurcations that are caused by the presence of singularities singular bifurcations, which is a more general consideration of possible scenarios, which may or may not involve equilibria directly. That is, we consider bifurcations caused by singularities including SIBs but not exclusively so.

To make a clear impression of the realm of all possible bifurcations, we focus on low-dimensional quasilinear DAEs of form (Equation 1) for which x ∈ ℝ or x ∈ ℝ². Our goal is to provide a list of all possible singular bifurcations of codimension 1 in such systems.

The paper is organized as follows. In Section 2 the basic notions for the study of one-dimensional quasilinear DAEs are given, all possible codimension-one bifurcations are studied and the behavior in higher-codimension bifurcations is described. In Section 3 the main notions are given for a two-dimensional case, for which the full list of possible codimension-one bifurcations is provided. Section 4 contains the rigorous derivation of the dynamical behavior near some local singular bifurcations from Section 3.

Consider a quasilinear DAE (Equation 1) for n = 1, it is given by (4)g(x,α)x.=f(x,α), (x,α)∈ℝ×ℝm, for smooth functions f : ℝ × ℝ^m → ℝ and g : ℝ × ℝ^m → ℝ. Note that in this case, the singular set is precisely the set of zeros of g given by (5)Σα={x∈ℝ:g(x,α)=0}, which under a regularity assumption on g, is composed of isolated points. We will call every such point a singularity. The set of zeros of f that are not zeros of g (6)Eα={x∈ℝ:f(x,α)=0,g(x,α)≠0}, will be called the equilibrium set. Under a similar regularity assumption, this set is also composed of isolated points. Each such point is an equilibrium of system (Equation 4).

Definition 2.1. A point x ∈ ℝ is called a singular equilibrium if it lies in the intersection of the singular set Σ_α given by Equation (5) with the zeros set of f(x, α) for some α ∈ ℝ^m.

Definition 2.2. For a fixed parameter value α = α₀, we will call the point x₀ a simple equilibrium if it is a simple zero of f and not a zero of g: f(x0,α0)=0,fx′(x0,α0)≠0,g(x0,α0)≠0. Analogously, we will call x₀ a simple singularity if it is a simple zero of g and not a zero of f: f(x0,α0)≠0,g(x0,α0)=0,gx′(x0,α0)≠0.

Definition 2.3. A simple singularity point x₀ is called incoming (outgoing), if there exists a small neighborhood U of x₀ such that for any initial condition x ∈ U the solution x(t) reaches x₀ in finite forward (backward) time.

Remark 2.4. The incoming or outgoing simple singularities are known in the DAE literature as the standard singular points which was introduced in [32]. They behave like an impasse point, where solutions are no longer defined being either attractive or repelling (cf. [32, 33]).

For a simple equilibrium point x₀ there exists a small neighborhood |x − x₀| < ε such that g(x, α₀) ≠ 0, and the system (Equation 4) can be rewritten as (7)x.=f(x,α0)g(x,α0):=f~(x,α0), where f~(x0,α0)=0. Thus, the stability type of the equilibrium (x₀, α₀) is completely determined by the sign of the derivative (8)λ:=∂xf~|(x0,α0)=g∂xf-f∂xgg2|(x0,α0)=∂xfg|(x0,α0)≠0, which is non-zero by the assumption that x₀ is a simple equilibrium for α = α₀. If λ < 0, then the equilibrium x₀ is stable; if λ > 0, then it is unstable (see Figure 1).

For a simple singularity x₀, there exists a small neighborhood in which f(x, α₀) ≠ 0 and the system (Equation 4) can be rewritten as (9)g~x.=1, for g~=g(x,α)f(x,α) with g~(x0,α0)=0. Thus, the following derivative (10)λ=g~x|(x0,α0)=gxf-gfxf2|(x0,α0)=gxf|(x0,α0)≠0, determines the type of the singularity at (x₀, α₀). More precisely, if λ > 0, then the simple singularity point x₀ is outgoing; if λ < 0, then it is incoming. See Figure 2, we put here and further below the double arrow to reflect the fact that the trajectory reaches the singularity in finite time, and the velocity grows to infinity.

It is clear from Definition 2.2, that simple equilibria and simple singularities persist under generic parametric perturbations. Indeed, the condition fx′(x0,α0)≠0 implies that by the Implicit Function Theorem, there exists locally a unique function x^*(α) with x*(α0)=x0, which fulfills the equation f(x, α) = 0. Moreover, this equilibrium maintains the same stability type, as the exponent λ in Equation (8) preserves its sign. In a similar way, one can deduct the corresponding property for a simple singularity.

Theorem 2.5. If in an open set U ⊂ ℝ system (Equation 4) possesses a finite set of equilibria and a finite set of singularities and all of them are simple, then the system is structurally stable in U.

Proof: Under sufficiently small perturbations, every simple equilibrium and every simple singularity stays in a small neighborhood of its initial position, remain simple, are distributed in the same order on the line and keep their stability types. Also, neither of these points reaches the boundary of U. The intervals bounded by these equilibria and singularities can be homeomorphically conjugated, with the direction of motion preserved.

However, when the conditions of Theorem 2.5 are violated, one may encounter bifurcations. There are three such possibilities of singular bifurcations¹:

A1. A non-simple equilibrium: f(x, α) = 0, fx′(x,α)=0;

These cases are not exclusive to each other. They can happen simultaneously, either at the same or different points, which increases the codimension of the problem. In the following, we assume that the bifurcation conditions occur at a₀ = 0 and x₀ = 0, which can be achieved by appropriate translation of coordinates and parameters.

We start with formulating the simplest possible cases, i.e. the cases of codimension-1.

The case A1.1 is completely analogous to the usual fold bifurcation in dynamical systems without singularities. Indeed, since g(x, α) ≠ 0 in a small neighborhood of (x₀, α₀) = (0, 0), we can rewrite the system (Equation 4) using the function f~ as defined in Equation (7), where f~(0,0)=f~x′(0,0)=0, f~″xx(0,0)=f″xxg|(0,0)≠0.

Proposition 2.6. Assume that for system (Equation 4) the conditions of case A1.1 are fulfilled for α = 0 at x = 0, and fα′(0,0)≠0. Then for all small α by an invertible change of coordinate and parameter, the system can be brought near the origin to the following normal form: (11)η.=β+sη2+O(η3), where s=signfxx″g|(0,0) (cf. Figure 3 left).

This result immediately follows from Kuznetsov [34], Theorem 3.2. Similarly, one can derive the normal form for the non-simple singularity in case A2.1.

Proposition 2.7. Assume that for the system (Equation 4) the conditions of case A2.1 are fulfilled for α = 0 at x = 0, and gα′(0,0)≠0. Then, for all small α by an invertible change of coordinate and parameter, the system can be brought near the origin to the following normal form: (12)(β+sη2+O(η3))η.=1, where s=signgxx″f|(0,0) (cf. Figure 3 right).

Proof: In some small neighborhood of the origin, we have f(x, α) ≠ 0 for all small α. Then, the system can be rewritten in the form (Equation 9) with g~(0,0)=g~x′(0,0)=0, g~″xx(0,0)=g″xxf|(0,0)≠0. We expand g~ in Taylor series in x as g~=g0(α)+g1(α)x+g2(α)x2+O(x3) with g₀(0) = g₁(0) = 0 and g₂(0) = a ≠ 0. Using a parameter-dependent coordinate shift of the form x = y + δ(α) with δ(0) = 0, we can rewrite g~ as: g~(y+δ(α),α)=(g0(α)+g1(α)δ(α)+O(α2))+(g1(α)+2g2(α)δ(α)+O(α2))y+(g2(α)+O(α))y2+O(y3), where the linear term can be neglected by an appropriate choice of δ. Indeed, the coefficient of the linear term vanishes at α = δ = 0, and its derivative with respect to δ at zero is given by 2a ≠ 0. By the Implicit Function Theorem, there exists a function δ(α)=δ1α+O(α2) with δ1=-g1,α′(0)2a. Thus, we have [g0,α′(0)α+O(α2)+a(α)y2+O(y3)]y.=1 with a(0) = a, which becomes (12) using the scaling y=1|a(α)|1/3η, β=1|a(α)|1/3(g0,α′(0)α+O(α2)). Moreover, as g0,α′(0)=gα′f(0,0)≠0, all the transformations are invertible.

The bifurcation occurs in the following way, if s = −1: for β = 0 there exists a locally unique non-simple singularity point, which disappears when β < 0 and is replaced by a pair of incoming and outgoing simple singularity points when β > 0 [see Figure 3 (right)].

The following Proposition states the normal form of the transcritical singularity bifurcation in case A3.0,0.

Proposition 2.8. If the system (Equation 4) satisfies the conditions of transcritical singularity, case A3.0.0, at (x, α) = (0, 0), and A:=(gx′fα′-fx′gα′)|(0,0)≠0, then there exists an invertible change of coordinate and parameter, which brings system (Equation 4) near the origin to the normal form (13)ηη.=β+sη+O(η2), where s=signfx′gx′|(0,0).

(cf. Figure 4).

Proof: As gx′(0,0)≠0, the implicit equation g(x, α) = 0 can be locally uniquely resolved with respect to x for small x and α. That is, there exists a smooth function x^*(α) such that g(x^*(α), α) ≡ 0 with x^*(0) = 0 and x*(α)=-gα′gx′|(0,0)α+O(α2). Consider a parameter-dependent shift of coordinate x = x^*(α) + y. Then, the left-hand side of Equation (4) is transformed as (14)g(x*(α)+y,α)=gx′(x*(α),α)y+O(y2)=(g1+O(α))y+O(y2)                           =(g1+O(α))y(1+O(y)), where g1=gx′(0,0)≠0. The right-hand side function f becomes (15)f(x*(α)+y,α)=f(x*(α),α)+fx′(x*(α),α)y+O(y2)                      =Ag1α+O(α2)+(f1+O(α))y+O(y2), where f1=fx′(0,0).

We choose the neighborhood small enough for term (1 + O(y)) in formula (Equation 14) to stay always positive. Then, we reparametrize time by formula dt/(1 + O(y)) = dτ, and also divide by a non-zero coefficient (g₁ + O(α)), system (Equation 4) is transformed as: (16)yy.=Ag12α+O(α2)+(f1g1+O(α))y+O(y2). It leads to Equation (13) by scaling y→|f1g1+O(α)|η and setting (17)β=Af12α+O(α2). Notice that as A ≠ 0, the new small parameter β diffeomorphically depends on α.

Example 2.9. The following systems demonstrate examples of cases A1.1, A2.1 and A3.0,0, respectively.

More precisely, consider the system (Equation 18) with f(x, α) = x² + α and g(x, α) = x + 1 at (x, α) = (0, 0). It satisfies the non-simple equilibrium conditions, case A1.1: f(0,0)=fx′(0,0)=0, f″xx(0,0)=2≠0, g(0,0)=1≠0, fα′(0,0)=1≠0. By Proposition 2.6, the normal form of this bifurcation is given by Equation (11) with s=signfxx″g|(0,0)=1.

For the system (Equation 19) at the non-simple singularity point (x, α) = (0, 0), one has g(0,0)=gx′(0,0)=0, g″xx(0,0)=2≠0, f(0,0)=1≠0, gα′(0,0)=1≠0. By Proposition 2.7, the normal form here is (12) with s=signgxx″f|(0,0)=1.

The system (Equation 20) at point (x, α) = (0, 0) satisfies the conditions of the transcritical singularity: f(0,0)=g(0,0)=0, fx′(0,0)=1≠0, gx′(0,0)=1≠0, |gx′gα′fx′fa′|(0,0)=1≠0. By Proposition 2.8, the normal form of this bifurcation is Equation (13) with s=signfx′gx′|(0,0)=1.

Besides the codimension-1 bifurcations listed by A1.1–A3.0,0, one can describe generic unfoldings of bifurcations of higher codimension in a similar way. These bifurcations admit invertible changes of coordinates (shifts), that bring them to the corresponding normal forms. We do not give a proof here, because it is straightforward: similar to cases A1.1 and A2.1 there always exists the parameter-dependent shift of coordinates x → y + δ(α) such that the term that vanishes at α = 0 and has the highest power in the Taylor expansion (y^m or yⁿ below), is eliminated for all small α. We distinguish the following bifurcations:

A1.m. A codimension-m equilibrium: f(0, 0) = 0, fx(i)(0,0)=0 for 1 ≤ i ≤ m and fx(m+1)(0,0)≠0, g(0, 0) ≠ 0. The normal form is: (21)y.=β0+β1y+β2y2+…+βm-1ym-1+sym+1+O(ym+2)

In the formulas above coefficient s is equal to either +1 or −1, and all α_i and β_i are small unfolding parameters.

Lemma 2.10. In one-dimensional system (Equation 4) the higher-order bifurcations occur in the way that under small perturbations the following dynamics is observed, depending on the case A1–A3.

Proof: First, we consider cases A1.m and A2.n. Take a set of integers {a₁, a₂, …, a_k} as described above in the respective case, and select small x₁ < x₂ < … < x_k (this means, that |x_i| < ε for some ε). Then, construct the following polynomial: (24)P(y)=(y-x1)a1(y-x2)a2…(y-xk)ak. It has k roots at coordinates x_i with multiplicity a_i each. If we open the parentheses in formula (24), the resulting polynomial will be a small perturbation of the order-(m + 1) polynomial standing in the right-hand side of normal form (21) or the order-(n + 1) polynomial in the left-hand side of normal form (22). The statement of Lemma on equilibria or singularities respectively, directly follows.

Now take case A3.m, n and the respective sets of integers {a₁, a₂, …, a_k} and {b₁, b₂, …, b_l} and coordinates: x₁ < x₂ < … < x_k and y₁ < y₂ < … < y_l. We construct two polynomials: (25)P(y)=(y-x1)a1(y-x2)a2…(y-xk)ak,Q(y)=(y-y1)b1(y-y2)b2…(y-yl)bl. Polynomial P(y) is the small perturbation of the right-hand side and polynomial Q(y) is the small perturbation of the left-hand side of normal form (23). At the same time, this system possesses equilibria, singularities and singular equilibria exactly as described in the Lemma.

Consider the two-dimensional quasilinear DAEs of form (Equation 1), where A is everywhere nonsingular except on the singular set Σ_α. The simplest possible form of such DAEs is given by (cf. [35]) (26){g(x,y,α)x.=f1(x,y,α)y.=f2(x,y,α), for (x, y) ∈ ℝ², α ∈ ℝ^m and g:ℝ² × ℝ^m → ℝ, f1,f2:ℝ2×ℝm→ℝ are smooth functions.

In this case, the singular set Σ={(x,y):g(x,y,α)=0} is the zero curve of g. This curve is the boundary of two domains: Σ+={(x,y):g(x,y,α)>0}, Σ-={(x,y):g(x,y,α)<0}. We assume that g is such that g_x has finitely many zeros on Σ. This means that ∇g is zero also at finitely many points of Σ. Every point of Σ with ∇g ≠ 0 belongs to the closure of both Σ⁺ and Σ⁻.

In order to describe the dynamics of such a two-dimensional system, we introduce the basic dynamical elements, such as special points and cycles.

Definition 3.1. A point (x, y) is called

Remark 3.2. The fold points given by Definition 3.1 are also referred to as non-standard algebraic singular points in DAE terminology (cf. [32, 33]). The singular equilibria are also called as non-standard in some literature, e.g., in [32], or standard (as extended in [33]) geometric singular points.

Making time transformation dτ:=dtg for g ≠ 0 in Equation (26), one obtains the desingularized system: (27){x.=f1(x,y,α)y.=f2(x,y,α)g(x,y,α). The correspondence between the systems is the following: every trajectory Γ of system (Equation 27) contains one or more trajectories of system (Equation 26). The intersection points of Γ with singularity curve g = 0 (if any) split Γ into connected pieces Γ₁, Γ₂, … that are trajectories of Equation (26) with the same direction of time as Γ for every Γi⊂Σ+ and with the opposite direction of time if Γi⊂Σ-. It is clear that both equilibria and singular equilibria of system (Equation 26), given by Definition 3.1, are equilibria of the ODE system (Equation 27) in the usual sense.

Remark 3.3. Transforming the time by formula dτ:=dt-g creates another desingularized system, in which the time flows in the opposite direction on every trajectory of Equation (27). The trajectories of system (Equation 26) follow the trajectories of this system in Σ⁻ and flow in the opposite direction in Σ⁺.

Consider a point M ∈ Σ and trajectory Γ ∋ M of ODE system (Equation 27), see [36] for details. The trajectory will intersect Σ transversely if it is not an equilibrium state and if its tangent vector (f1,f2·g)⊤|M=(f1,0)⊤ is not orthogonal to ∇g=(gx,gy)⊤, i.e. f₁ · g_x ≠ 0. Locally, Γ is split by M into two components, Γ₁ and Γ₂ such that M = Γ(0), Γ₁ ⊂ Γ(t) for t < 0 and Γ₂ ⊂ Γ(t) for t > 0. If f₁ · g_x > 0, then Γ crosses Σ from Σ⁻ to Σ⁺, Γ2⊂Σ+ is the trajectory of Equation (26), and Γ1⊂Σ- is the trajectory of Equation (26) with time reversal. In this case we call point M outgoing. When f₁ · g_x < 0, the direction of time is preserved on Γ1⊂Σ+ and reversed on Γ2⊂Σ-, and point M is called incoming. Inequalities f₁ · g_x > 0 and f₁ · g_x < 0 are open conditions, thus curve Σ consists of incoming Σ^inc and outgoing Σ^out zones, separated by points where f₁ · g_x = 0, those are either singular equilibria f₁ = 0 or fold points g_x = 0.

Definition 3.4. A limit cycle of system (Equation 27) is called a limit cycle of system (Equation 26), if it has no intersections with the singularity curve Σ. Otherwise, it is called a folded limit cycle.

The limit cycle is a periodic orbit of system (Equation 26). A folded limit cycle consists of more than one orbit of system (Equation 26).

In this section the character of those local bifurcations, that do not have analogs in regular ODEs is studied in detail. These bifurcations are L3: Singular saddle-node of type II, L4: Cubic fold and L5: Singular equilibrium-fold from subsection 3.2.2.

We introduce the following notations: (36)Δ1=det∂(f1,f2)∂(x,y)|(0,0,0), Δ2=det∂(f1,g)∂(x,y)|(0,0,0), Δ3=det∂(g,gx)∂(y,α)|(0,0,0), Δ4=∂(f1,f2,g)∂(x,y,α)|(0,0,0), Δ5=∂(f1,g,gx)∂(x,y,α)|(0,0,0).