Transient Responses to Relaxation Oscillations in Multivibrators

The multivibrator is an electronic circuit that has three oscillation states: astable, monostable, and bistable; these circuits typically contain opamps. These states are often modeled using hybrid systems, which contain characteristics of both continuous and discrete time. While an ideal opamp possesses both continuous and discrete characteristics, actual opamps exhibit continuous properties, which necessitate in-depth modeling. The relaxation oscillations produced by the multivibrator, characterized by periodic, rapid state changes, are typically modeled by considering slow–fast dynamical systems. In these systems, the phenomenon whereby the amplitude of the signal changes rapidly is referred to as a “canard explosion”. By considering this phenomenon, it is possible to understand the process of relaxation oscillations in the multivibrator. In this work, we model the multivibrator by considering a slow-fast dynamical system and observe canard explosions through numerical experiments. This study indicates that the oscillatory changes in the multivibrator are continuous, which explains the onset of relaxation oscillations. Additionally, circuit experiments are conducted using affordable opamps; in this experimental work, canard explosions are observed.


I. INTRODUCTION
A multivibrator [1] is a type of electronic circuit typically implemented using opamps.These circuits are often used as timers or switches.Multivibrators, as the name suggests, can generate multiple types of oscillations.The oscillation states of the system vary depending on the circuit configuration and are classified into three types: 1. Astable: This variety of multivibrator continuously produces oscillations and is commonly used as an oscillator circuit.2. Monostable: This system produces oscillations once and then stops.3. Bistable: Bistable multivibrators have two stable states, and the oscillation state that is active is based on the initial conditions or external inputs.
These multivibrators are often modeled as hybrid systems.A hybrid system possesses characteristics of The associate editor coordinating the review of this manuscript and approving it for publication was Ming Xu .both continuous-time and discrete-time dynamical systems.In multivibrators that are assumed to contain an ideal opamp, the state not only undergoes continuous changes but also experiences discrete changes due to switching.When such circuits are modeled as a hybrid system, one can neglect the transient responses of the oscillation states, making it relatively easy to derive return maps.
The Astable mode of the multivibrator is often used as a square wave oscillator.In related research on dynamical systems, electronic fireflies [2] and their synchronization phenomena [3] have been analyzed using square wave oscillators.In these studies, the square wave oscillators are examined as hybrid systems using ideal operational amplifiers, and precise return maps have been obtained.However, the opamps and operational amplifiers (opamps) used in the realization of these circuits do not possess ideal characteristics.Thus, the system must be modeled as a continuous system when we consider actual circuit VOLUME 12, 2024VOLUME 12, 2023 The Authors.This work is licensed under a Creative Commons Attribution 4.0 License.
For more information, see https://creativecommons.org/licenses/by/4.0/systems.Even sophisticated opamps and opamps produce outputs with slight delays [1].The transient responses of the opamps within multivibrators should not be overlooked.
Monostable and bistable multivibrators are designed such that their equilibrium state is also an equilibrium point of the system (rather than an oscillation state being an equilibrium point), making the analysis of state changes relatively straightforward.However, since astable systems continuously produce oscillations, the delay in the opamp induces different oscillating states.It is particularly apparent in astable multivibrators that transitions from non-oscillating to oscillating states occur when the circuit parameters change.At these transitions, despite the actual device having continuous characteristics, the amplitude of the circuit output is observed to change ''discontinuously'' from zero.One might intuitively assume that this amplitude explosion is continuous [4].However, the transient response of the stable state to such parameter changes has not been investigated previously.
The square waves that are produced by multivibrators are often referred to as relaxation oscillations [5].Relaxation oscillations involve a rapid state change within a certain cycle, followed by a period in which that state is maintained.Thus, ''relaxation'' implies that both slow and fast changes are involved in a transition between states.Such oscillations can be explained by considering a slow-fast dynamical system [4].These are systems composed of two continuous-time dynamical systems that operate on different timescales; these systems permit both slow-and fast-state characteristics to be investigated.
A notable phenomenon observed in slow-fast dynamics is the canard [6].Canards occur only in a very limited range of parameters immediately prior to an amplitude explosion [4], [7].A canard is a solution that changes its amplitude drastically in response to a small parameter change; this kind of amplitude explosion is referred to as a canard explosion.We also note that the term ''Canard'' is derived from the French word for ''duck'', and it refers to the characteristic ''duck-like headed shape'' of the trajectory with large amplitudes.When it can be shown that the dynamics of a multivibrator contains a canard, it is possible to demonstrate the continuity of the amplitude change that occurs during the relaxation oscillation.
This work models a multivibrator as a slow-fast dynamical system and numerically investigates the canard explosion that occurs due to parameters changing during the relaxation oscillation.It is found that the transition to the relaxation oscillation of the multivibrator is continuous.We provide an explanation for this observation based on the bifurcation theory of dynamical systems [8], [9].Additionally, based on the results obtained via the numerical experiments presented here, we conduct circuit experiments.Economical opamps are used to observe canards easily, demonstrating the possibility of observing this complex phenomenon without the use of expensive equipment.In the experiments, both canards and canard explosions are observed.

II. SLOW-FAST DYNAMICAL SYSTEMS
A slow-fast dynamical system [4] can be represented by a system of ordinary differential equations of the form: where and 0 < ϵ ≪ 1. Dividing both sides of the second equation in (1) by ϵ, dy/dt becomes larger than dx/dt.Therefore, in this work, x is referred to as the slow variable and y is called the fast variable.By setting τ = t/ϵ, we can also obtain the equivalent form of (1): In this work, we refer to the dynamics produced by Eqs.( 1) and ( 2) as the slow-timescale and the fast-timescale dynamics, respectively Let us consider the case of the singular limit ϵ → 0. In this limit, (1) becomes a differential-algebraic equation, which is given by, and ( 2) becomes a layer equation, We refer to (3) and ( 4) as the reduced problem.The set defined by the second equation in (3), is referred to as the critical manifold.The flow described in (3) can be considered to be determined by dx/dt subject to the condition that it is bounded on C 0 .Fig. 1 shows an schematic illustration example of a typical flow within slow-fast dynamical systems for m = n = 1.The figure is the case of van der Pol equation which is a most classical [6] slow-fast dynamical system.We consider f to be cubic function which is a type of y 3 and g to be a linear function.The points p − and p + in C 0,s = {p ∈ C 0 : ∂g/∂y(p, 0) is not invertible} represent points where the uniqueness of the flow is lost; these are called fold points and satisfy the following expressions: In the singular limit ϵ → 0, C 0 can be divided into subsets according to C 0 = C + 0,a ∪ C 0,r ∪ C − 0,a ∪ C 0,s .Of these subsets, C 0,a = C + 0,a ∪ C − 0,a represents the attractive part of the critical manifold and C 0,r represents the repelling part of C 0 .If an equilibrium point lies on C 0,a , the equilibrium point is completely stable.When the equilibrium point is situated on C 0,r , the equilibrium point is completely unstable and there exists a stable periodic solution.It can be seen that the periodic solution evolves along C ± 0,a and jumps to C ∓ 0,a when the orbit reaches a folding point.This behavior is called a relaxation oscillation.The dynamics of this transition can be considered to be a hybrid system based on (3) and ( 4), with the reaching of the folding point being considered as the event.
One of the notable phenomena in slow-fast dynamical systems is the canard solution [6].The canard solution occurs when the equilibrium point on C 0,r transitions to a relaxation oscillation via a Hopf bifurcation.In such a system, the canard is observed only within a very small region in the parameter space, and a amplitude of a stable periodic solution increases explosively as a result of a small change in the parameters.This behavior is referred to as a canard explosion.Fig. 2 shows a schematic one-parameter bifurcation diagram around a Hopf bifurcation.This figure is drawn based on the illustration used in Fig. 8.3 of [4].In this figure, the parameter A denotes the amplitude of the attractor, A = max y − min y, where max y, min y are the maximal and minimal value of the limit cycle.As depicted in Fig. 2, canard explosions can be classified into two types: (a) those associated with a supercritical Hopf bifurcation and (b) those associated with a subcritical Hopf bifurcation.Solid lines and dotted lines represent stable and unstable limit cycles/equilibria, respectively.In the case of Fig. 2(b), the tangent bifurcation is observed at the parameter where stable and unstable periodic solutions adhere together.This indicates that the system is bistable when the parameter is between the tangent bifurcation and the Hopf bifurcation.
The van der Pol oscillator represents a typical slow-fast dynamical system [6].In the original van der Pol equation [10], a vacuum tube amplifier is used, resulting in a critical manifold of sigmoid shape.However, we consider a model using a cubic function for simplicity: Fig. 3 shows an example of limit cycles and canard solutions in a van der Pol equation for ϵ ̸ = 0 (see ( 6)).Fig. 3(a1-a4) show the phase portraits of the system for ϵ = 1 and q = 1.05, 0.987, 0.9863, and 0.5, respectively.This system does not exhibit slow-fast characteristics.The amplitude of the trajectory changes continuously as the parameter changes.Fig. 3(b1-b4) are the phase portraits of the system for ϵ = 0.1 and q = 1.05, 0.987, 0.9863, and 0.5, respectively.It is noted that slow-fast characteristics can be observed in this system for ϵ = 0.1.In all the figures, regardless of the initial state, the trajectories converge to the orbits shown in the figures for those parameters.We note that the change in the value of q used to obtain Figs.3(b2) and ( b3) is very small, and this small change in q leads to large changes in the behavior of the system.This rapid amplitude change represents a canard explosion.We classify the small-amplitude cycle immediately before the observed canard explosion (Figs.3(b2)) as a ''canard without head'' and the large-amplitude cycle immediately after the explosion (Figs.3(b3)) as a ''canard with head.''It can be observed that even a small change in the parameters describing the system can lead to a rapid increase in the amplitude.In the case of van der Pol oscillator, Analytical methods for calculating the Canard explosion parameter are provided [11].Fig. 4 shows the amplitude changes that are induced by parameter variations in a van der Pol oscillator.In Fig. 4, the grey lines q ± indicate the value of q at which the equilibrium points coincide with the fold points p ± .From the figure, it can be seen that when the equilibrium point is on C 0,r , a Hopf bifurcation occurs and periodic solutions emerge.Conversely, when the equilibrium point is on C 0,a , the amplitude is zero, indicating that no periodic solutions occur.In systems that correspond to the behavior shown in Fig. 4(a), the amplitude of the stable limit cycle changes smoothly.However, in the slow-fast dynamical system characterized by Fig. 4(b), the amplitude increases abruptly.Typically, as the value of ϵ decreases, the amplitude explosion becomes increasingly steep; this makes the rise in amplitude appear to be discontinuous.
Another typical system in which canards can be observed is the FitzHugh-Nagumo model [12].The FitzHugh-Nagumo FIGURE 3. Stable equilibrium points and limit cycles in a van der Pol oscillator.(a1-a4) depict the case of ϵ = 1, and (b1-b4) show the case of ϵ = 0.1.The trajectories obtained for the various parameters and an example of a canard are shown.We refer to (b2) as a canard without head and (b3) as a canard with head.FIGURE 4. Amplitude changes in a van der Pol oscillator.The system without slow-fast dynamics, shown in subfigure (a), exhibits a smooth amplitude change, whereas in the case of the system that exhibits slow-fast dynamics (subfigure (b)), the amplitude increases abruptly.The grey line indicates the parameter values at which the equilibrium points coincide with the fold points p ± .model describes the electrical activity of neurons, and the rapid changes resembling spike responses can be attributed to the slow-fast dynamics of the system.In the coupled FitzHugh-Nagumo model, canards are observed, and the canard explosion is analytically determined [13].Canards can also be observed in discrete-time spiking neuron dynamics [14] and self-replicating systems [15].We also note the existence of an interesting canard phenomenon in aircraft trajectories reported in Ref. [16].

III. MULTIVIBRATOR
In this section, we obtain a circuit of a multivibrator as a slow-fast dynamical system.First, we construct a multivibrator as a hybrid system using an ideal operational amplifier and explain how square wave oscillations are generated.Next, we demonstrate that by using the dynamic characteristics of the operational amplifier, the system is constituted as a slow-fast dynamical system.Thereafter, we describe the relationship between the modes of the multivibrator and its equilibrium points, providing the prerequisite knowledge necessary for subsequent numerical calculations.We consider the circuit shown in Fig. 5.
Fig. 6 shows the output characteristics of a single-supply opamp.The output of an ideal opamp takes the form of a step function with an increase at v d = 0, where v d is the input voltage difference, given by 474 VOLUME 12, 2024 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.By using v d , the output of the opamp is expressed as v o = a(v d ).In the case of an ideal opamp, a(v d ) becomes a step function.The actual output characteristics of an opamp take the shape of a sigmoid curve, but first, we consider the dynamics when using an ideal opamp.
According to Kirchhoff's law, the relationship between the input voltage, v p and v o , can be described as follows: where, we replace coefficients of (8) as: Then, we have the ideal opamp output as: For the RC circuit at the bottom of the opamp, we have: From ( 7) to (10) and the output characteristic with step function v o = a(v d ), we can obtain the system as the following differential algebraic equation: Fig. 7 shows the dynamics of the multivibrator system with opamp.For sake of simplicity, we consider the case of R C → ∞.In this limit, the slope of f (x, y) is equal to 1. Here, we consider the case where v d = 0, that is, the point at which a(v d ) is discontinuous.In the case of v o = v p − v C = E, we have v C = (β + γ )E.Similarly, in the case of v o = 0, we have v C = γ E. When the system is considered as a hybrid system, these values represent the points where the system triggers an event.The dynamics will be explained using points a, b, c, and d in the figure.Suppose the initial value is given at point a.In this case, the state v C changes according to the differential equation (11).Next, when v C reaches b, a(v d ) becomes 0 and makes a discontinuous jump to c.This corresponds to an event trigger.Subsequently, following the differential equation, v C decreases and reaches d.When the system reaches d, a(v d ) = E and it returns to point a.This results in the relaxation oscillation of the multivibrator.Fig. 8 shows the time-domain response in Fig. 7.In the figure ,  points a, b, c, and d   In practice, opamps have an output characterized by a steep increase (but not infinitely steep) around v d = 0; this output has a form similar to that of a sigmoid function.The static output characteristic of the opamp can be described as Considering the dynamic output characteristic as a first-order lag system [17], we obtain: where τ 0 is a time constant and a parameter that causes the circuit to behave as a slow-fast dynamical system.In this work, we approximate the output characteristic of the opamp using the hyperbolic tangent function according to: where α represents the gain.Note that an ideal opamp is characterized by τ 0 = 0 in (12) and α = ∞ in (13).
From ( 8), ( 10), ( 12) and ( 13), we obtain the following system of second-order differential equations: where v c and v o are replaced with x and y, respectively; these substitutions are performed to obtain a notation consistent with (1).Additionally, we re-scaled the time constant by RCt =: t and we set ϵ = τ 0 /RC to match the form to (1).
The parameter ϵ includes the parameters R and C, thus for example, changing R in equation ( 14) will also change ϵ.
In this study, we will fix these values of R and C. Note that the changing RC can adjust the slow-fast dynamics easily.Moreover, by adjusting the value of R C , the slope of the linear equation on the right-hand side of the first equation of ( 14) can be varied, allowing the position of the equilibrium points to be easily manipulated.Unless otherwise noted, in this work, we set and Fig. 9 shows a classification of the locations in which the equilibrium points can be generated.In the case of an equilibrium point being generated at C 0,a , as shown in 1 ⃝ and 3 ⃝ in Fig. 9, the equilibrium point is completely stable.On the other hand, when an equilibrium point is generated on C 0,r , as shown in 2 ⃝ in Fig. 9, it may become an unstable equilibrium point via a Hopf bifurcation; the classification of this equilibrium point depends on the value of ϵ.This corresponds to the cause in which relaxation oscillations are observed.In other words, the characteristics of the multivibrator can be adjusted by changing the position of the equilibrium point.
In our proposed model, by varying the value of R C , the slope of f (x, y) is changed.This results in the multivibrator depicted in Fig. 5 exhibiting two modes: one mode remains at a stable equilibrium point, whereas the other mode exhibits rectangular oscillations.
It is interesting to consider when the equilibrium point is near one of the two singular points, p ± ; we show an example illustrating this situation in Fig. 10.As indicated by one of the dashed lines in the figure, is it possible for the trajectory to exhibit small amplitudes?The relaxation oscillation of the multivibrator starts drawing suddenly large amplitudes as shown by the solid line.During the process of the equilibrium point moving from C − 0,a to C 0,r , a square-wave oscillation suddenly emerges.This implies that despite the continuous variation in the parameters defining the system, the amplitude changes discontinuously.This ''transient response due to parameter variation'' can be explained by considering the canard explosion in slow-fast dynamical systems.

IV. NUMERICAL ANALYSIS
In slow-fast dynamical systems, due to the disparity in the timescales, ordinary numerical integration methods often suffer from a loss of accuracy or even generate fake chaotic trajectories (trajectories that do not actually exist) [18].Analytical methods using singular perturbation theory [19] are available, but in this study, a classical numerical method of the dynamical systems e.g. a shooting method and a numerical continuation method are used, and we discuss the existence of canard explosions from the perspective of bifurcation theory.It is thus necessary to use appropriate numerical integration methods and utilize methods such as multi-precision arithmetic in order to ensure sufficient accuracy.Here, we use the Runge-Kutta-Fehlberg method to perform the necessary numerical integration.
In this section, we compute numerically the canard explosion points [20] and canard solutions.On the faster timescale, we consider the following planar system [7]: where f and g are C ∞ -class function, λ ∈ R is a parameter, and 0 < ϵ ≪ 1.We assume that C 0 is locally parabolic and the minimum is coincident with the origin, (0, 0).We call the origin a fold point, and this point satisfies: for λ ̸ = 0. We can obtain the slow flow on C 0 by differentiating x = ϕ(y) with respect to t = τ ϵ: where the function x = ϕ(y) for ϕ : U → R, U is sufficiently small neighborhood of y = 0.The slow flow is singular at the origin for λ ̸ = 0 since dϕ/dt(0) = 0 and f (0, 0, λ, 0) ̸ = 0. We assume a non-degenerate canard point, which is an equilibrium point located on the origin, for λ = 0.A canard point satisfies the following additional conditions: f (0, 0, 0, 0) = 0, ∂f ∂x (0, 0, 0, 0) ̸ = 0, ∂f ∂λ (0, 0, 0, 0) ̸ = 0. ( This gives us a well-defined slow flow on C 0 for λ ̸ = 0. Near a non-degenerate canard point, we have a normal form [21]: where We obtain the Hopf bifurcation parameter λ H and the canard explosion parameter λ c as [7]: where K H and K c is are real numbers defined by h 1 -h 6 .Refer to the reference [7] for a detailed definition of K H and K c .It can be seen that the equilibrium is stable for λ < λ H and unstable for λ > λ H .The type of Hopf bifurcations can be identified by considering the sign of K c ; supercritical bifurcations exist for K c < 0, and K c > 0 indicates a subcritical bifurcation.We see that another expression relating λ H and λ c can be obtained: We can obtain the Hopf bifurcation parameter, λ H , via conventional numerical methods.Thus, given K c , the canard explosion parameter λ c can be obtained.
When the system is written in a form with strong nonlinearity, it is difficult to apply the method of transformation to the normal form described here.Therefore, to avoid equation transformations, Kuehn developed a method to numerically calculate λ c using the first Lyapunov coefficient [20].The first Lyapunov coefficient, l 1 , is equal to K c scaled by a constant, meaning we can obtain λ c by computing l 1 numerically.However, there are several definitions of the first Lyapunov method [8], [22].This is due to the background in traditional dynamical systems theory, where the sign of the first Lyapunov coefficient is important, and the actual value of the coefficient does not need to be considered.Therefore, the scaling factor ρ will change depending on the type of first Lyapunov coefficient used.In this work, we use Kuznetsov's convention and notate it l Ku 1 [8].
We assume that the equilibrium point (x * , y * ) of ( 2) is under Hopf bifurcation and is translated to coincide with the origin with the coordinate change z = (x − x * , y − y * ) ⊤ , and we thus obtain: dz dt = M z + F(z) (24) with F(z) = O(||z|| 2 ) and M ∈ R (m+n)×(m+n) .This form (24) is the linearization around the equilibrium point of (2), M is the Jacobian matrix around the equilibrium point.Taking a Taylor series expansion of the nonlinear term F, we have, where the multilinear functions B and C are defined as, where B(u, v) and C(u, v, w) are symmetric multilinear vector functions of u, v, w ∈ R (m+n) , F i denotes the i-th element of the function F. In the case of a planar system, we thus obtain a simple form of l Ku 1 : where ω 0 is given by the eigenvalues of the matrix M , λ 1,2 = ±iω 0 , g 20 = p⊤ B(q, q), g 11 = p⊤ B(q, q), and g 21 = p⊤ C(q, q, q).ℜ takes the real part of complex number.p, q ∈ C (m+n) are eigenvectors of λ 1 and the transpose M ⊤ , respectively.These are chosen such that they satisfy p⊤ q = 1.
Note that M has the pure imaginary eigenvalues 0 + iω 0 since we consider the equilibrium point which undergoes Hopf bifurcation.
The first Lyapunov coefficient, l Ku 1 , has the following property [4], [20]: where ρ is the positive scaling factor.We then obtain an expression for λ c , where ρ = 1/ ρ.In (30), we can obtain the scaling factor ρ by calculating λ H and λ c .λ H can be obtain by numerical bifurcation analysis which we show later at (31).Also, we can obtain approximated values of λ c with numerical continuation method.We show the specific scaling factor ρ in the numerical analysis later.In this paper we treat R C as the role of λ.
The Hopf bifurcation parameter, λ H , is obtained by solving the following conditions to obtain the parameters of an equilibrium point and λ H , (x * , y * , λ H ), numerically: where J denotes the Jacobian matrix of (f , g) ⊤ with respect to (x, y) ⊤ , ⊙ denotes the bialternate product [23], and I is the (m + n) × (m + n) identity matrix.Here, we use Newton's method to solve the Hopf bifurcation condition.
We show the actual procedure to obtain the canard explosion parameter λ c below.
1) Obtain the Hopf bifurcation parameter λ H .We use the shooting method with (31).2) Calculate the first Lyapunov coefficient l 1 at the Hopf bifurcation.3) Compute the maximal canard parameter λ c with the (30).Note that the scaling factor ρ is obtained in advance by numerically as we described before.To obtain the scaling factor ρ, we use ϵ = 0.001.In this paper we use the approximated parameter λ c = 104130.3954at ϵ = 0.001 which is obtained from continuation of periodic orbits using Runge-Kutta-Fehlberg method, the Hopf bifurcation parameter λ H = 105543.3018with ω 0 = 0.157901 which is obtained from the shooting method with (31).The first Lyapunov coefficient at the Hopf bifurcation is l Ku 1 = 0.87939173.Then we have ρ = 203231491.0351and able to calculate the canard explosion set with (30).This scale ρ changes with the parameter selected as λ.In this case, we have chosen the parameter  R C , which takes on large values, resulting in a large scale as well.
In slow-fast dynamical systems, it is difficult to compute λ c ; this calculation typically involves numerical integration, and it is difficult due to the precision requirements of the integrator because of the slow-fast characteristics [18].Kuehn's first Lyapunov coefficient method [20] offers the advantage of involving only algebraic operations and permits the computation of λ c without utilizing numerical integration.Fig. 11 shows the (R C , 1/ϵ) bifurcation diagram.The set represented by the curve labelled H corresponds to the Hopf bifurcation, while C ± represent the sets of canard explosion points near the singular points, p ± .The gray lines in the figure represent the parameter values at which the equilibrium and singular points coincide.The locations of p ± are obtained by solving the following condition for (x * , y * , λ 0 ) by Newton's method: where λ 0 is the parameter which the equilibrium point coincide to fold point.In this case, it is We obtain the equilibrium point (x * , y * ) which coincides to a fold point and the parameter λ 0 at the same time by solving the objective function.We show the actual values: R − C = 106.0540876[k] and its location is (x * , y * ) = (x 0 , y 0 ) = (0.7282, 0.3748), R + C = 11.1570095[k] and its location is (x 0 , y 0 ) = (4.2717,0.4287), where (x 0 , y 0 ) shows a fold point.In the parameter region within the Hopf bifurcation through C ± , unstable equilibria and stable periodic solutions emerge.The stable periodic solutions are characterized by small amplitudes immediately after the Hopf bifurcation, but the canard explosion at C ± induces a rapid increase in their amplitude.Further changes in the parameter values lead to relaxation oscillations.
Fig. 12 shows the one-parameter bifurcation diagram of the amplitudes of periodic solutions; this figure corresponds to 1/ϵ = 300 in Fig. 11.The amplitude A is the difference of the maximal and minimal values of limit cycles as same as Fig. 2. The gray lines in Fig. 12, labelled R ± C , represent the parameter that a equilibrium point coincides to a fold point, and the magenta lines C ± represent the canard explosion parameter, as is the case in Fig. 11.The amplitude can be seen to increase rapidly when the Hopf bifurcation occurs.It appears in Fig. 12 that the precision of C + is poor, but, as mentioned above, the precision improves as ϵ decreases.Indeed, small values of ϵ lead to sharp increases in the amplitude.The relaxation oscillation observed in the multivibrator can then be attributed to canard explosions that are present for sufficiently small values of ϵ.
The classification of canards as ''canards with head'' and ''canards without head'' can be determined by whether the limit cycle contains points with zero curvature or not [24].Consider the planar system given in (1).Trajectories can be obtained by eliminating time, t, from the equation, Differentiating this expression with respect to y gives, The trajectories with zero curvature thus satisfy, By plotting the set that satisfies (35) on the x-y plane, it is possible to determine whether the limit cycle has points with zero curvature.Fig. 13 shows periodic orbits near the canard explosion.Fig. 13(a) represents a ''canard without head'', and Fig. 13(c) shows a ''canard with head''.The parameters for Figs.13(a The curves labelled I ϵ in Fig. 13 represent the set of solutions where the curvature of the trajectory is 0. It can be seen that the headless canard does not intersect with I ϵ , whereas the headed canard does.Fig. 13(b) represents a case close to the canard explosion point.Here, we observe a fake chaotic trajectory, which occurs due to the slow-fast dynamics.This is because the accuracy of the numerical integration decreases.The multivibrator model proposed in this study is a two-dimensional autonomous system, and such trajectories are not permissible.These fake chaotic trajectories will be confirmed in the subsequent circuit experiments.

V. CIRCUIT IMPLEMENTATION
In the above numerical simulations, we observed that the multivibrator considered here undergoes a transition from a stable state to a relaxation oscillation, and canards are observed during this process.The region in the parameter space in which canards can be observed is small, but canards have been observed [25] in multiple nonlinear electronic circuits [26].To validate the findings of the numerical work presented here, we implemented the multivibrator in a circuit and demonstrated the occurrence of canard explosions experimentally.Here, we aim to demonstrate the occurrence of canards in the real circuit response of the multivibrator and capture the topological changes due to the canard explosion, a characteristic of slow-fast dynamical systems.Since our goal is not to precisely replicate the canards as shown in numerical calculations, we do not verify errors in components such as resistance elements.
Due to the large value of α in (13), the actual components of the opamp induce a sharp Z-shape in the curve of C 0 , and the amplitude explosion after the Hopf bifurcation is very pronounced.Furthermore, we note that the time-delay characteristic corresponding to ϵ is determined by the slew rate of the opamp; this value is sufficiently small.In the numerical calculations presented in the previous section, R G was set to be smaller than R E and R F to reduce the severity of the Z-shape and to make the change in the slope of C 0 near p ± more gradual.This makes it possible to observe canards over a relatively wide range of parameter values.
Here, the circuit is implemented according to Fig. 5, using an Analog Devices OP177 opamp, which has a relatively low slew-rate; this means the system is more susceptible to canards than it would be using high slewrate opamps.This low-cost opamp has a gradual output characteristic, which makes it suitable for confirming canards.We note that it would be very challenging to observe canard explosions in systems containing highperformance opamps.However, even in the case of idealized opamps, numerical simulations have predicted the existence of canards within a very limited range of parameters.
Fig. 14 shows the circuit response of an experimentally realized multivibrator.It is interesting to observe the variations in the response of the circuit for different values of ϵ.However, since we cannot directly control the time constants of the opamp, we adjust C to modify the slow-fast characteristics of the system.The rows of the figure show the variation in the observed trajectories for a fixed value of C as a result of changes in R C .As R C changes, the system transitions from being characterized by a canard without head via a canard explosion to being characterized by a canard with head (relaxation oscillations).The columns in Fig. 14 show the variations in the trajectories that occur as a result of changes in the value of C for a fixed value of R C .As C increases, the movement along the slow-variable direction becomes slower, resulting in more pronounced slow-fast characteristics, while the fast movement along the fast-variable direction becomes more emphasized.Both the headless and headed canards shown in the previous section are visible in (a2) and (b2) due to the slight parameter variations induced by small external noise.The multivibrator is a planar system, thus, such trajectories are not possible in the absence of external noise.

VI. CONCLUSION
In this work, we constructed a multivibrator as a slow-fast dynamical system, which exhibited responses typical of slow-fast systems.The proposed multivibrator allows for easy adjustment of the slow-fast characteristics by modifying the circuit components.Furthermore, the position of the equilibrium point of the system can be easily changed.The Hopf bifurcation set was obtained via conventional numerical computations.The set of canard explosion points was obtained using a method based on the first Lyapunov coefficient which only requires non-complicated numerical computation.Both the numerical simulations and circuit experiments undertaken here demonstrate the existence of canard explosions in the system.As future work, the series of numerical methods presented in this paper will be applied to other slow-fast dynamical systems.This includes applications to higher-dimensional systems where the fast dynamics involve two or more.In these higher-dimensional systems, analytical solutions are challenging, underscoring the importance of numerical computations.

FIGURE 1 .
FIGURE 1. Schematic illustration which represents a critical manifold and folding points.This illustration is based on the type of van der Pol equation.

FIGURE 2 .
FIGURE 2. Schematics of one-parameter bifurcation diagrams.Canard explosions, i.e., the explosive amplitude changes, are shown around the Hopf bifurcation parameter λ H . (a) Represents the case of a supercritical Hopf bifurcation and (b) shows the case of a subcritical Hopf bifurcation.This figure is drawn based on a Fig. 8.3 from [4].

FIGURE 5 .
FIGURE 5. Circuit diagram of the multivibrator considered here.

FIGURE 6 .
FIGURE 6.Output characteristics of an opamp.
Fig.7shows the dynamics of the multivibrator system with opamp.For sake of simplicity, we consider the case of R C → ∞.In this limit, the slope of f (x, y) is equal to 1. Here, we consider the case where v d = 0, that is, the point at which a(v d ) is discontinuous.In the case of v o = v p − v C = E, we have v C = (β + γ )E.Similarly, in the case of v o = 0, we have v C = γ E. When the system is considered as a hybrid system, these values represent the points where the system triggers an event.The dynamics will be explained using points a, b, c, and d in the figure.Suppose the initial value is given at point a.In this case, the state v C changes according to the differential equation(11).Next, when v C reaches b, a(v d ) becomes 0 and makes a discontinuous jump to c.This corresponds to an event trigger.Subsequently, following the differential equation, v C decreases and reaches d.When the system reaches d, a(v d ) = E and it returns to point a.This results in the relaxation oscillation of the multivibrator.Fig.8shows the time-domain response in Fig.7.In the figure, points a, b, c, and d corresponding to Fig.7are marked.It can be observed that the mode transitions in a way that switches for both v C and v o .

FIGURE 9 .
FIGURE 9. Classification of the positions of the equilibrium and singular points, p ± .

FIGURE 10 .
FIGURE 10.An schematic illustration of the situation when the equilibrium point is near one of the two singular points.

FIGURE 13 .
FIGURE 13.Examples of canards in a multivibrator: (a) a canard without head, (b) fake chaos near the canard explosion point, (c) a canard with head.In (b), due to the strong slow-fast dynamics, the accuracy of the numerical integration decreases, leading to the observation of fake chaos.The data shown in this figure is obtained for ϵ = 0.001.
) and (c) are R C = 104.131[k] and R C = 104.130[k], respectively, with ϵ = 0.001.It can be seen that the amplitude increases significantly even though the slope of f (x, y) = 0 changes only very slightly.

FIGURE 14 .
FIGURE 14.Canards observed in an experimental circuit.We capture the data using an Agilent DSO1024A oscilloscope.The sampling rate used to obtain the data presented here was 12.5 kHz.8192 data points are plotted in each subfigure.