Higher-Order Hamiltonian for Circuits with (α,β) Elements

Abstract

The paper studies the construction of the Hamiltonian for circuits built from the (α,β) elements of Chua’s periodic table. It starts from the Lagrange function, whose existence is limited to Σ-circuits, i.e., circuits built exclusively from elements located on a common Σ-diagonal of the table. We show that the Hamiltonian can also be constructed via the generalized Tellegen’s theorem. According to the ideas of predictive modeling, the resulting Hamiltonian is made up exclusively of the constitutive relations of the elements in the circuit. Within the frame of Ostrogradsky’s formalism, the simulation scheme of Σ-circuits is designed and examined with the example of a nonlinear Pais–Uhlenbeck oscillator.

1. Introduction

The dynamics of complex nonlinear systems attracts the interest of researchers in various branches of science [1,2,3,4]. One of the approaches used for studying system dynamics, predictive modeling, considers the system an interconnection of one-port devices, the α,β elements, also known as Higher Order Elements (HOEs) [5]. Each HOE guarantees that its terminal quantities are coupled in all circumstances via a constitutive relation. In electrical engineering, the terminal quantities are the generalized voltage v^(α) and current i^(β), where the integers α and β denote the orders of the derivative/integral with respect to time (for positive/negative integers). Well-known (α,β) elements include resistors (0,0), capacitors (0,−1), and inductors (−1,0), and also the promising memristors (−1,−1), memcapacitors (−1,−2), and meminductors (−2,−1), frequency-dependent resistors (FDNR, FDNC), and others. HOEs are usually represented as points with integer (α,β) coordinates in Chua’s periodic table (see Figure 1). The constitutive relation of the element is frequently written in one of the two following forms:

$$v^{\left(\alpha \right)}=f\left(i^{\left(\beta \right)}\right),\mathrm{or}{i}^{\left(\beta \right)}=g\left(v^{\left(\alpha \right)}\right)$$

(1)

where the function f() or g() models the current or voltage representation of the element.

The use of HOEs considerably increases the possibilities of analyzing and synthesizing circuits with complex dynamic behavior. Thanks to predictive modeling, many phenomena can be better understood. For example, the hypothetical element called a memristor [6] helped to explain the peculiar hysteresis behavior observed when measuring various nanomaterials. When Hewlett-Packard reported the first TiO₂-based working memristive device in 2008 [7], they reminded the researchers of the original work [6] that explains this behavior as a characteristic feature of the memristor—a fingerprint. Ongoing research has demonstrated that a pinched hysteresis loop is a regular phenomenon, which takes effect for an arbitrary (α,β) element in the space of the time derivatives of its terminal quantities [8].

The interest in systems with high-order dynamics has led to an extension of Lagrange’s and Hamilton’s formalisms to systems with α,β elements. This is because the original Lagrange’s formalism starts from the fundamental schematic of classical mechanics, i.e., from the point of view of predictive modeling based on the R (0,0), L (−1,0), and C (0,−1) elements. Moreover, the system of equations of motion can be generated from one scalar function, the Lagrangian, but only for conservative systems. Recent works [9,10,11] have extended Lagrange’s formalism to circuits containing memristors, memcapacitors, and meminductors. The following work [12] defines the potential functions for general (α,β) elements from which Lagrangians and dissipative functions are drawn, and rules for generating the equations of motions are established. The potential functions of a general (α,β) element are defined for current or voltage representation as integrals:

$$S_{\alpha ,\beta }={\displaystyle \int v^{\left(\alpha \right)}d}{i}^{\left(\beta \right)}\mathrm{or}{\widehat{S}}_{\alpha ,\beta }={\displaystyle \int i^{\left(\beta \right)}d}{v}^{\left(\alpha \right)},$$

(2)

which represent a natural generalization of the potential and kinetic energy and the dissipative function of classical mechanics, or the energy of capacitors, inductors, and the content of resistors in classical electrical engineering. Two different representations of potential functions correspond to the concept of functions and co-functions introduced by Millar [13] and Cherry [14] in 1951.

The classical Lagrangian is of the first order, so it is a function of generalized coordinates and velocities (first-order time derivatives of the coordinates); the resulting equations of motion are of the second order. Higher-order Lagrangians were introduced by Ostrogradsky in [17] as follows:

$$L=L\left(\mathit{x},{\mathit{x}}^{\left(1\right)},\dots ,{\mathit{x}}^{\left(m\right)},t\right)$$

(3)

The Lagrangian is a function of the 1xn vector of the generalized coordinates x = [x₁ .. x_n]^T and their derivatives up to the integer-order m. The system represented by Lagrangian (3) is governed by n equations of motion:

$$\frac{\partial L}{\partial \mathit{x}}-\frac{d}{dt}\left(\frac{\partial L}{\partial {\mathit{x}}^{\left(1\right)}}\right)+\frac{d^2}{d{t}^2}\left(\frac{\partial L}{\partial {\mathit{x}}^{\left(2\right)}}\right)-\dots +{\left(-1\right)}^m\frac{d^m}{d{t}^m}\left(\frac{\partial L}{\partial {\mathit{x}}^{\left(m\right)}}\right)=\mathbf{0}.$$

(4)

Equation (4) is a direct consequence of the fact that the system trajectory is extremal to the action

$$A={\displaystyle \underset{t_1}{\overset{t_2}{\int }}L\left(\mathit{x},{\mathit{x}}^{\left(1\right)},\dots ,{\mathit{x}}^{\left(m\right)},t\right)dt},$$

(5)

which is the essence of Hamilton’s variational principle [18]. The higher-order Lagrangian (3) was used for the first time to describe circuits built from HOEs in [16]. This work demonstrated that the Lagrangian, which has the ability to generate equations of motion, can only be drawn for the Σ-circuits that contain only elements from the common Σ-diagonal of Chua’s table, where the sum of the indices α and β is preserved; thus, Σ = α + β. The Lagrangian then has the form

$$L={\displaystyle \sum _{{\epsilon }_i}{\left(-1\right)}^i{S}_i}\mathrm{or}\widehat{L}={\displaystyle \sum _{{\epsilon }_i}{\left(-1\right)}^i{\widehat{S}}_i}.$$

(6)

The summation is done through all ε_i elements of the given Σ-diagonal, where i is the position of the element on the diagonal (see Figure 2). The signs in front of the state functions in the sums of (6) are governed by the type (even or odd) of the positions of the elements on the diagonal [16]. The difference between kinetic and potential energy, transferring between the inertial and the accumulating elements, which are the common form of the Lagrangian used in classical mechanics, are merely a special case of (6). The physical dimension of the Lagrangian is given by the number Σ of the diagonal: [Volt⋅Amper⋅sec^−Σ]. This is the energy [Volt⋅Amper⋅sec] in LC and the power [Volt⋅Amper] in resistive circuits. An example of the Lagrangian of the Pais–Uhlenbeck oscillator [15], consisting of three resistive elements of R, FDNR, and FDNC types, is given in [16].

The transition from the higher-order Lagrangian (3) to the higher-order Hamiltonian is given by the generalized Legendre transformation published by Ostrogradsky in [17]:

$$\mathcal{H}\left(\mathit{x},{\mathit{x}}^{\left(1\right)},‥,{\mathit{x}}^{\left(2m-1\right)}\right)={}^1\mathit{p}^T\cdot {\mathit{x}}^{\left(1\right)}+\dots +{}^m\mathit{p}^T\cdot {\mathit{x}}^{\left(m\right)}-L\left(\mathit{x},{\mathit{x}}^{\left(1\right)},\dots ,{\mathit{x}}^{\left(m\right)}\right)$$

(7)

where ¹p to ^mp are vectors of the generalized momenta

$${}^j\mathit{p}={\displaystyle \sum _{k=j}^m{\left(-1\right)}^{k-j}}{\left(\frac{\partial L}{\partial {\mathit{x}}^{\left(k\right)}}\right)}^{\left(k-j\right)}.$$

(8)

Consider the notation

$${}^1\mathit{q}=\mathit{x},{}^2\mathit{q}={\mathit{x}}^{\left(1\right)},\dots ,{}^m\mathit{q}={\mathit{x}}^{\left(m-1\right)}.$$

(9)

Assume that the Lagrangian is regular, i.e.,

$$\det \left(\frac{{\partial }^{2}L}{\partial x_i^{\left(m\right)}\partial x_j^{\left(m\right)}}\right)\ne 0,$$

(10)

for i,j = 1,...,n [19]. In Equation (8), for j = m reads as

$${}^m\mathit{p}=\frac{\partial L}{\partial {\mathit{x}}^{\left(m\right)}}.$$

(11)

Then, the unambiguous relation,

$${\mathit{x}}^{\left(m\right)}=\mathit{\chi }\left({}^1\mathit{q},\dots ,{}^m\mathit{q},{}^m\mathit{p}\right),$$

(12)

can be derived from (11), which enables a transition to new coordinates ¹p,.., ^mp and ¹q,..., ^mq and allows us to draw the canonic equations,

$${}^j\dot{\mathit{q}}=\frac{\partial H}{\partial {}^j\mathit{p}},{}^j\dot{\mathit{p}}=-\frac{\partial H}{\partial {}^j\mathit{q}}.$$

(13)

The Hamiltonian H is then a function of the new variables ⁱq, ⁱp, i = 1, ... , m:

$$H\left({}^1\mathit{q},\dots ,{}^m\mathit{q},{}^1\mathit{p},\dots ,{}^m\mathit{p}\right)={}^1\mathit{p}^T\cdot {}^2\mathit{q}+\dots +{}^{m-1}\mathit{p}^T\cdot {}^m\mathit{q}+{}^m\mathit{p}^T\cdot \mathit{\chi }\left({}^1\mathit{q},\dots ,{}^m\mathit{q},{}^m\mathit{p}\right)-\mathcal{L}\left({}^1\mathit{q},\dots ,{}^m\mathit{q},{}^m\mathit{p}\right)$$

(14)

where

$$\mathcal{L}\left({}^1\mathit{q},\dots ,{}^m\mathit{q},{}^m\mathit{p}\right)=L\left({}^1\mathit{q},\dots ,{}^m\mathit{q},\mathit{\chi }\left({}^1\mathit{q},\dots ,{}^m\mathit{q},{}^m\mathit{p}\right)\right).$$

(15)

The procedure of drawing an alternative Hamiltonian for cases when condition (10) of the regularity of the Lagrangian is not fulfilled is described, for example, in [19]. Irrespective of whether the Lagrangian is or is not regular, the Hamiltonian can be used within Lagrange’s formalism, i.e., in the sense of generalized energy (7).

This formalism, specific to higher-order Hamiltonians of type (14), is denoted Ostrogradsky’s formalism [19]. None of the hitherto published works deal with assembling the Hamiltonian of systems consisting of elements from Chua’s table. The benefits of using the Hamiltonian have been proven: the system dynamics is given by a set of canonic first-order differential equations, which are the standard starting point for solving complex tasks associated with terms such as the Lyapunov exponents, zero divergence, conservative phase volumes [20], etc. Hamiltonians are an ideal tool when searching for symmetries and associated conservative quantities [21]. Today, they are well defined also for systems that have hitherto been considered non-conservative [22].

The objective of this work is therefore to find the Hamiltonian of a Σ-system consisting of arbitrary HOEs within the frame of both Lagrange’s formalism (as generalized energy) and Ostrogradsky’s formalism. The conventional procedure for introducing Ostrogradsky’s formalism into physical theory involves searching for the form of the Lagrangian and generalized momenta, and subsequently obtaining from them a specific form of the Hamiltonian ((7) or (14)). The procedure used in this work will be physically more objective, starting from the generalized Tellegen’s theorem. In addition to its clarity, the reason for selecting this approach is pragmatic: Tellegen’s theorem is extremely general and applicable to circuits of arbitrary topologies with any elements (linear, nonlinear, time-invariant, and time-varying), as well as with arbitrary elements from Chua’s table. Moreover, this theorem can also be formulated via the terminal quantities of individual elements. It is useful for circuits containing HOEs, since the sole fundamental characteristic of the element—the constitutive relation—is a link only between these terminal quantities. Then, the resulting Hamiltonian will naturally be comprised of the constitutive relations of all HOEs in the circuit.

This paper has the following structure. Section 2 introduces the generalized form of Tellegen’s theorem applicable to circuits with HOEs and thus to circuits with the elements defined via the constitutive relations (1). Section 3 is devoted to drawing a specific form of the Hamiltonian for Σ-circuits in the frame of Lagrange’s formalism and also of Ostrogradsky’s formalism. Section 4 describes a modeling technique inspired by Ostrogradsky’s formalism. In the last section, these new pieces of knowledge are applied to a specific topology of the Pais–Uhlenbeck oscillator consisting of HOEs.

2. Tellegen‘s Theorem for Circuits with Higher-Order Elements

Consider a circuit comprised of arbitrary one-port elements ε with port voltages v_ε and currents i_ε. Then, the generalized Tellegen’s theorem [23] holds for such a circuit:

$${\displaystyle \sum _{\epsilon }{v}_{\epsilon }^{\left(\alpha \right)}{i}_{\epsilon }^{\left(\beta \right)}}={\left({\mathit{v}}_{\epsilon }^{\left(\alpha \right)}\right)}^T\cdot {\mathit{i}}_{\epsilon }^{\left(\beta \right)}=0$$

(16)

where α and β are arbitrary integers, ${\mathit{v}}_{\epsilon }^{\left(\alpha \right)}$ and ${\mathit{i}}_{\epsilon }^{\left(\beta \right)}$ are 1xb vectors of generalized voltages and currents of the element, and b is the number of elements in the circuit. The classical form of the theorem, representing the case α = β = 0, means that the sum of instantaneous powers delivered to all elements in the circuit is zero. The generalized theorem (16) replaces the instantaneous power [VA] by a quantity [VA]^{α + β}. This is because the integration and differentiation with respect to time belong to Kirchhoff’s operators [19], which do not affect the validity of Kirchhoff’s laws and the theorems derived from them. Note that the source of the voltage v^(α) or current i^(β) can be substituted by an (α,β) element with a constant constitutive relation (1) of a f() or g() type, so theorem (16) also holds for circuits with voltage and current sources.

Consider a circuit comprised of general HOEs. The element ε_h = (α_max,β_min) is the hidden element of the current representation of the circuit [24]. Let us introduce new variables u = v^(αmax) and x = i^(βmin). Then, Equation (16) can be rewritten with the aid of Equation (1) in the form

$${\displaystyle \sum _{\epsilon }{f}_{\epsilon }^{\left(\alpha +{\alpha }_{\max }-{\alpha }_{\epsilon }\right)}{x}_{\epsilon }^{\left(\beta +{\beta }_{\epsilon }-{\beta }_{\min }\right)}}=0$$

(17)

where f_ε and x_ε are the constitutive relation (1) and the x variable of the element ε. Since Equation (16) must hold for arbitrary integers α and β, it also holds for α = 0 and β = 1 + β_min − β_ε. Substituting these values in formula (17), which is identically equal to zero, and integrating it with respect to time, the resulting quantity must be constant in all circumstances:

$$\mathcal{H}={\displaystyle \int {\displaystyle \sum _{\epsilon }{f}_{\epsilon }^{\left({\mathsf{\Delta }}_{\epsilon }\alpha \right)}{\dot{x}}_{\epsilon }}dt}=const$$

(18)

where Δ_εα = α_max − α_ε is the distance of the ε element from the hidden element in the α direction (see Figure 2a).

A similar procedure leads to the dual form of Equation (18) for voltage representation of the circuit:

$${\mathcal{H}}^{*}={\displaystyle \int {\displaystyle \sum _{\epsilon }{g}_{\epsilon }^{\left({\mathsf{\Delta }}_{\epsilon }\beta \right)}{\dot{u}}_{\epsilon }}dt}=cons{t}^{*}$$

(19)

where Δ_εβ = β_ε − β_min is the distance of the ε element from the hidden element ε_h = (α_min,β_max) in the β direction (see Figure 2b).

It will be shown in the following section that the quantity (18) or (19) is, for the case of Σ-circuits, a Hamiltonian of type (7) and that it becomes a Hamiltonian of type (14) after transitioning to coordinates that correspond to Ostrogradsky’s formalism.

3. From Tellegen‘s Theorem to the Hamiltonian

Let all the elements of the circuit be located on a common Σ-diagonal according to Figure 2. Denote such a circuit as the Σ-circuit. All types of elements will be specified as ε₀ to ε_m, where m is the distance of the ε_m element from the hidden element in the α or β direction (since the element is located on the Σ-diagonal, these distances are the same, thus Δ_εα = Δ_εβ =Δ_ε). The method of indexing depends on the position of the hidden element, and this position depends on the choice between the current, f(), and the voltage, g(), representation (1) of the element. The relation ε∈ε_i, i=0,...,m will thus specify that the given ε element is only of the ε_i type.

Let us focus on the proof that the formula (18) for the current representation of elements (see Figure 2a) is a Hamiltonian. The proof for the formula (19), which corresponds to the voltage representation, is analogous, so, for the sake of brevity, it will not be given below.

The sum (18) for ε∈ε₀ is equal to

$${\displaystyle \sum _{\epsilon \in {\epsilon }_0}{\displaystyle \int f_{\epsilon }^{\left(0\right)}{\dot{x}}_{\epsilon }dt}}={\displaystyle \sum _{\epsilon \in {\epsilon }_0}{\displaystyle \int f_{\epsilon }^{\left(0\right)}d{x}_{\epsilon }}}={\displaystyle \sum _{\epsilon \in {\epsilon }_0}{S}_{\epsilon }}.$$

(20)

For ε∈ε_m, m > 0, repeated integration by parts leads to the following result:

$${\displaystyle \sum _{\epsilon \in {\epsilon }_m}{\displaystyle \int f_{\epsilon }^{\left(m\right)}{\dot{x}}_{\epsilon }dt}}={\displaystyle \sum _{\epsilon \in {\epsilon }_m}\left({\displaystyle \sum _{j=0}^{m-1}{\left(-1\right)}^j{f}_{\epsilon }^{\left(m-1-j\right)}{x}_{\epsilon }^{\left(j+1\right)}}\right)}+{\displaystyle \sum _{\epsilon \in {\epsilon }_m}{\left(-1\right)}^m{S}_{\epsilon }}.$$

(21)

Summing for all the circuit elements and rearranging the terms yield the conservative quantity (18) in a compact form:

$$\mathcal{H}={\displaystyle \sum _{j=1}^m{\left(-1\right)}^{j-1}\left({\displaystyle \sum _{i=j}^m{\displaystyle \sum _{\epsilon \in {\epsilon }_i}{x}_{\epsilon }^{\left(j\right)}{f}_{\epsilon }^{\left(i-j\right)}}}\right)}+{\displaystyle \sum _{i=0}^m{\displaystyle \sum _{\epsilon \in {\epsilon }_i}{\left(-1\right)}^i{S}_{\epsilon }}}.$$

(22)

The generalized currents x_ε can be represented as linear combinations of generalized loop currents x. Utilizing incidence matrices, the addend of the inner summation of the first term in (22) will assume the form

$${\displaystyle \sum _{k=1}^n{a}_{\epsilon k}{x}_k^{\left(j\right)}{f}_{\epsilon }^{\left(i-j\right)}}={\displaystyle \sum _{k=1}^n{x}_k^{\left(j\right)}{}^{k}f_{\epsilon }^{\left(i-j\right)}}$$

(23)

where ^kf_ε means either ±f_ε or 0 depending on whether the element ε is or is not a part of the k-th loop, and, possibly, depending on what its orientation is with regard to this loop. Substituting (23) into (22), expanding the outer series according to index j, and subsequently rearranging the summation yield

$$\mathcal{H}={\displaystyle \sum _{k=1}^n\left(x_k^{\left(1\right)}\underset{{}^{1}p_k}{\underbrace{{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{\epsilon \in {\epsilon }_i}\left({}^{k}f_{\epsilon }^{\left(i-1\right)}\right)}}}}+x_k^{\left(2\right)}\underset{{}^{2}p_k}{\underbrace{{\displaystyle \sum _{i=2}^m{\displaystyle \sum _{\epsilon \in {\epsilon }_i}\left(-{}^{k}f_{\epsilon }^{\left(i-2\right)}\right)}}}}+\dots +x_k^{\left(m\right)}\underset{{}^{m}p_k}{\underbrace{{\displaystyle \sum _{\epsilon \in {\epsilon }_m}({\left(-1\right)}^{m-1}{}^{k}f_{\epsilon }^{\left(0\right)})}}}\right)}-\underset{L}{\underbrace{{\displaystyle \sum _{i=0}^m{\displaystyle \sum _{\epsilon \in {\epsilon }_i}{\left(-1\right)}^i{S}_{\epsilon }}}}}$$

(24)

Since it holds that

$$\frac{\partial L}{\partial x_k}={\left(-1\right)}^{k+1}{}^{k}f_{\epsilon },$$

(25)

comparing (24) and (7) reveals that the function $\mathscr{H}$ from Equation (24) is the generalized energy (7) of the current representation of a circuit with HOEs, and the generalized momenta are

$${}^{j}p_k={\left(-1\right)}^{j-1}{\displaystyle \sum _{i=j}^m{\displaystyle \sum _{\epsilon \in {\epsilon }_i}\left({}^{k}f_{\epsilon }^{\left(i-j\right)}\right)}}.$$

(26)

The one-dimensional case for n = 1 signifies that all the elements are in series, with the common generalized current x. The series connection of the elements of the same type can, therefore, be regarded as one element of the same type, with the constitutive relation given as the sum of constitutive relations of the individual elements. Equation (26) can, therefore, be rewritten in the simplified form

$${}^{j}p={\left(-1\right)}^{j-1}{\displaystyle \sum _{i=j}^m{f}_i^{\left(i-j\right)}\left(x^{\left(i\right)}\right)}$$

(27)

where f_i() denotes the constitutive relation of an element of the ε_i type. The generalized momenta are then given as summations of the time derivatives of the constitutive relations

$$\begin{gathered} {}^{1}p=f_1^{\left(0\right)}+f_2^{\left(1\right)}+‥+f_m^{\left(m-1\right)} \\ {}^{1}p=f_1^{\left(0\right)}+f_2^{\left(1\right)}+‥+f_m^{\left(m-1\right)} \\ ⋮ \\ {}^{m}p={\left(-1\right)}^{m-1}{f}_m^{\left(0\right)}. \end{gathered}$$

(28)

This one-dimensional case clearly illustrates the transition from the formulation of the Hamiltonian (7) to the formulation (14) by changing the coordinates (x⁽⁰⁾,..., x^(m)) of Lagrange’s formalism to new coordinates (¹p,...,^mp, ¹q,...,^mq) of Ostrogradsky’s formalism. The key relation χ() from (12) can be obtained via a simple inversion of the constitutive relation of the element of the ε_m type:

$${}^{m}p={\left(-1\right)}^{m-1}{f}_m\left(x^{\left(m\right)}\right)\Rightarrow x^{\left(m\right)}=g_m\left({\left(-1\right)}^{m-1}{}^{m}p\right)$$

(29)

where g_m() is the constitutive relation of the voltage representation of the element of the ε_m type.

4. Modeling of the Σ-Circuits

Ostrogradsky’s formalism leads to a system of canonic first-order equations (13). This fact facilitates the modeling process. It follows from (28) that the momenta of a one-dimensional system are

$$\begin{gathered} {}^1\dot{p}=-f_0\left({}^{1}q\right) \\ {}^2\dot{p}=-{}^{1}p+f_1\left({}^{2}q\right) \\ {}^3\dot{p}={}^{2}p+f_2\left({}^{3}q\right) \\ ⋮ \\ {}^m\dot{p}={\left(-1\right)}^{m-1}{}^{m-1}p+f_{m-1}\left({}^{m}q\right). \end{gathered}$$

(30)

Similarly, the generalized coordinates are

$$\begin{gathered} {}^1\dot{q}={}^{2}q \\ {}^2\dot{q}={}^{3}q \\ ⋮ \\ {}^{m-1}\dot{q}={}^{m}q \\ {}^m\dot{q}=g_m\left({\left(-1\right)}^m{}^{m}p\right). \end{gathered}$$

(31)

The systems of Equations (30) and (31) lead to the elegant programming diagram in Figure 3, consisting of 2xm integrators for computing m generalized momenta and m generalized coordinates, as well as m+1 function blocks for modeling the constitutive relations of the individual elements.

5. Application: Non-Linear Pais–Uhlenbeck Oscillator

The usefulness of the Hamiltonian for the analysis of HOE circuits can be demonstrated by the example of the Pais–Uhlenbeck (PU) oscillator [15]. The PU oscillator is frequently employed to test new physical theories. It generates a signal with two harmonic components whose amplitudes and initial phases are given by the initial conditions. Both frequencies are coupled via the formula [25]

$${\omega }^{\pm }=\omega \sqrt{\frac{1\mp \sqrt{1-4ϵ}}{2ϵ}}$$

(32)

where ϵ is a real number, 0 ≤ ϵ < ¼. For ϵ → 0, the PU oscillator changes to a classical sinusoidal oscillator with an oscillation frequency of ω.

The PU oscillator is governed by the differential equation

$$\frac{ϵ}{{\omega }^2}\stackrel{⃜}{x}+\ddot{x}+{\omega }^{2}x=0.$$

(33)

In general, the signal x generated by the oscillator is non-periodical. Periodicity is achieved in special cases when the ratio of both frequencies is a rational number. The corresponding values of the constant ϵ are as follows:

$$ϵ=\frac{1}{4}\left(1-{\left(\frac{1-k^2}{1+k^2}\right)}^2\right),k=\frac{{\omega }^{+}}{{\omega }^{-}}=\frac{r}{s}<1$$

(34)

where r, s are integers. Only in these cases will the state-space trajectories be closed curves.

It is shown in [16] that the PU oscillator, modeled by Equation (33), can be implemented as a series or parallel Σ-circuit consisting of an arbitrary triad of HOEs, which are immediate neighbors on an arbitrary Σ-diagonal. One possible combination is shown in the inset in Figure 4a, namely, the series connection of the FDNC, R, and FDNR elements with the constitutive relations f₀(), f₁(), and f₂(). In general, these relations can be nonlinear. The indices 0 to 2 correspond to the distance of the element from the hidden element in the α or β direction. The hidden element is the FDNC. The generalized current x is governed by the differential equation

$${\ddot{f}}_2\left(\ddot{x}\right)+{\dot{f}}_1\left(\dot{x}\right)+f_0\left(x\right)=0$$

(35)

If all three elements have linear constitutive relations,

$$f_0\left(x\right)={\omega }^{2}x,f_1\left(\dot{x}\right)=\dot{x},f_2\left(\ddot{x}\right)=\frac{ϵ}{{\omega }^2}\ddot{x},$$

(36)

then the generalized current will be modeled via the linear differential Equation (33).

The subsequent analysis is made for linear FDNC and R and nonlinear FDNR elements with the following constitutive relations:

$$f_0\left(x\right)=k_{0}x,f_1\left(\dot{x}\right)=k_1\dot{x},f_2\left(\ddot{x}\right)=k_2\arctan \ddot{x}.$$

(37)

The Lagrangian of the circuit is

$$L\left(x,\dot{x},\ddot{x}\right)=-\frac{1}{2}{k}_0{x}^2+\frac{1}{2}{k}_1{\dot{x}}^2-\frac{1}{2}{k}_2\left(2\ddot{x}\arctan \ddot{x}-\ln \left(1+{\ddot{x}}^2\right)\right).$$

(38)

The equation of motion, generated by the Lagrangian according to (4), takes the form

$$k_{0}x+k_1\ddot{x}+k_2\frac{\stackrel{⃜}{x}\left(1+\ddot{x}\right)-{\stackrel{⃛}{x}}^2}{{\left(1+\ddot{x}\right)}^2}=0.$$

(39)

According to (28), the generalized momenta are

$$\begin{gathered} {}^{1}p=k_1\dot{x}+k_2\frac{\stackrel{⃛}{x}}{1+{\ddot{x}}^2} \\ {}^{2}p=-k_2\arctan \ddot{x}. \end{gathered}$$

(40)

For Lagrange’s formalism, the Hamiltonian (7) will be in the form

$$\mathcal{H}\left(x,\dot{x},\ddot{x},\stackrel{⃛}{x}\right)=\frac{1}{2}{k}_0{x}^2+\frac{1}{2}{k}_1{\dot{x}}^2+k_2\frac{\dot{x}\stackrel{⃛}{x}}{1+{\ddot{x}}^2}-\frac{1}{2}{k}_2\ln \left(1+{\ddot{x}}^2\right).$$

(41)

In conformity with (9), we introduce the new variables

$${}^{1}q=x,{}^{2}q=\dot{x}.$$

(42)

Utilizing Equation (40) to formulate the second derivative of x via the momentum ²p and slightly rearranging the equation will lead to the Hamiltonian

$$H\left({}^{1}q,{}^{2}q,{}^{1}p,{}^{2}p\right)={}^{1}p\cdot {}^{2}q+\frac{1}{2}{k}_0\cdot {}^{1}q^2-\frac{1}{2}{k}_1\cdot {}^{2}q^2+k_2\ln \left(\cos \left(\frac{{}^{2}p}{k_2}\right)\right).$$

(43)

The programming diagram for the nonlinear Pais–Uhlenbeck oscillator is shown in Figure 5.

The results of a computer simulation are shown in Figure 6. Independent variables as arguments of the generalized energy $\mathscr{H}$ (x, x⁽¹⁾, x⁽²⁾, x⁽³⁾) (41) and the Hamiltonian H (¹q, ²q, ¹p, ²p) (43) are graphically differentiated.

The parameters of the oscillator and initial conditions are selected such that the oscillation performs at the fundamental frequencies of 1 Hz and 7 Hz. The phase trajectories are, therefore, closed curves. For the sake of completeness, the value of ϵ does not correspond to the frequency ratio 1:7, as provided by formula (34) because (34) holds for the linear PU oscillator, while in the experiment, a nonlinear FDNR was used.

Figure 7 displays the waveforms of the Lagrangian and the Hamiltonian. While the Lagrangian is time-varying, the Hamiltonian, as a conservative quantity, holds its initial value, which is given by the initial conditions (see the legend in Figure 6). According to Equation (41) or (43), this value is 1.775 mW, which conforms with the simulation. Due to time discretization, the numerical integration causes a seeming loss of circuit conservatism [27]. The inset in Figure 7 shows negligible deflections of the computed Hamiltonian from its theoretical constant value.

The phase trajectories in Figure 8 represent the contours of the generalized energy (41): The condition$\mathscr{H}$(x, x⁽¹⁾, x⁽²⁾, x⁽³⁾) = const holds during the motion along these trajectories.

The Hamiltonian (43) remains constant during motion. The phase trajectories from Figure 9, therefore, represent the contours of the Hamiltonian H: The condition H(¹q, ²q, ¹p, ²p) = const holds during the motion along these trajectories.

Another possible implementation of the PU oscillator is shown in Figure 4b. The oscillator is built from the memristor MR and the linear (0,−2) and (1,−3) elements. Utilizing the MOVE transformation [24] in Chua’s table and the duality principle, the two last mentioned elements can be replaced by their linear versions of FDNC and FDPC (Frequency Dependent Positive Conductance). The coordinates of the hidden element (FDPC) indicate that the generalized current x is a threefold time integral of the current, or i⁽⁻³⁾. Equation (35) now expresses the KV⁽¹⁾L, where f₀, f₁, and f₂ are the constitutive relations of the FDPC, FDNR, and MR elements.

Consider linear FDPC and FDNR elements in series with the HP (Hewlett-Packard) memristor [7], with the nonlinear dopant drift being modeled via the Joglekar window function with the parameter p = 1 [28]. The constitutive relations will be in the form [29]

$$f_0\left(x\right)=k_{0}x,f_1\left(\dot{x}\right)=k_1\dot{x},f_2\left(\ddot{x}\right)=R_{off}\ddot{x}-\frac{\mathsf{\Delta }R}{4k}\ln \left(\frac{A\exp \left(4k\ddot{x}\right)+1}{B}\right)$$

(44)

where k₀, k₁ are positive constants, ΔR = R_off−R_on is the difference between the maximum and the minimum memristance, A = (R_off−R_ini)/(R_ini−R_on), B = ΔR/(R_ini−R_on), R_ini is the initial memristance, and k is a technological parameter [28]. The Lagrangian of the oscillator is generated after the integration of the constitutive relations and their appropriate summation:

$$L\left(x,\dot{x},\ddot{x}\right)=-\frac{1}{2}{k}_0{x}^2+\frac{1}{2}{k}_1{\dot{x}}^2+\frac{1}{2}{R}_{off}{\ddot{x}}^2+\frac{\mathsf{\Delta }R}{4k}\ddot{x}\ln B+\frac{\mathsf{\Delta }R}{{\left(4k\right)}^2}L{i}_2\left(-A\exp \left(4k\ddot{x}\right)\right)$$

(45)

where Li₂() denotes the dilogarithm [30], a special form of a polylogarithm, or Jonquière’s function.

The generalized momenta (28) are now

$$\begin{gathered} {}^{1}p=k_1\dot{x}+\left(R_{on}+\frac{\mathsf{\Delta }R}{A\exp \left(4k\ddot{x}\right)+1}\right)\stackrel{⃛}{x}. \\ {}^{2}p=-R_{off}\ddot{x}+\frac{\mathsf{\Delta }R}{4k}\ln \frac{A\exp \left(4k\ddot{x}\right)+1}{B} \end{gathered}$$

(46)

For Lagrange’s formalism, the Hamiltonian (7) is

$$\begin{gathered} \mathcal{H}\left(x,\dot{x},\ddot{x},\stackrel{⃛}{x}\right)=\frac{1}{2}{k}_0{x}^2+\frac{1}{2}{k}_1{\dot{x}}^2-\frac{1}{2}{R}_{off}{\ddot{x}}^2+\frac{\mathsf{\Delta }R}{4k}\ddot{x}\ln \left(A\exp \left(4k\ddot{x}\right)+1\right)+ \\ +\frac{\mathsf{\Delta }R}{{\left(4k\right)}^2}L{i}_2\left(-A\exp \left(4k\ddot{x}\right)\right)+\left(R_{on}+\frac{\mathsf{\Delta }R}{A\exp \left(4k\ddot{x}\right)+1}\right)\dot{x}\stackrel{⃛}{x}. \end{gathered}$$

(47)

Applying Ostrogradsky’s formalism, the Hamiltonian acquires the form

$$\begin{gathered} H\left({}^{1}q,{}^{2}q,{}^{1}p,{}^{2}p\right)={}^{1}p\cdot {}^{2}q+\frac{1}{2}{k}_0\cdot {}^{1}q{}^2-\frac{1}{2}{k}_1\cdot {}^{2}q^2+ \\ +\frac{\mathsf{\Delta }R}{{\left(4k\right)}^2}L{i}_2\left(-A\exp \left(4k{g}_2\left({}^{2}p\right)\right)\right)+g_2\left({}^{2}p\right)\left({}^{2}p+\frac{\mathsf{\Delta }R}{4k}\ln B+\frac{1}{2}{R}_{off}{g}_2\left({}^{2}p\right)\right) \end{gathered}$$

(48)

where g₂() is, according to (31), the inverse of the constitutive relation f₂() of the memristor.

The modeling diagram is the same as in Figure 5, but the functional block in the feedback will now correspond to the function g₂(). The inverse of f₂(), according to (44), cannot be derived in terms of standard functions. However, the inversion can be done numerically during the simulation run in SPICE. It follows from (44) that

$$\ddot{x}=h(f_2,\ddot{x})=\frac{f_2}{R_{off}}+\frac{\mathsf{\Delta }R}{4k{R}_{off}}\ln \left(\frac{A\exp \left(4k\ddot{x}\right)+1}{B}\right).$$

(49)

Based on the known value f₂, the simulation program computes $\ddot{x}$ from Equation (49) via iterations according to the SPICE code in

Table 1. The SPICE code for the numerical inversion of the constitutive relation $f_2(\ddot{x})$ in (44). For example, for f₂ = 120.9 (see the last line of the code), SPICE computes the voltage at node xdd, V(xdd) = 4 mV, which is equal to the value of the charge $\ddot{x}$. For simplicity, definitions of the parameters of function h are omitted. The function logexp prevents exponential overflow in (49).

.func logexp(x,A) {if(abs(x)<50,ln(1+A*exp(x)),ln(A)+x)}

.func h(f2,xdd)={xdd/Roff+delta/(4*k*Roff)*(−ln(B)+logexp(4*k*f2,A))}

Exdd xdd 0 value={h(v(xdd),120.9)}

.

The simulation results are summarized in Figure 10.

The parameters of the elements and simulation options are again set such that the generated waveforms are periodical. The phase trajectories are, therefore, closed curves.

Figure 11 displays the waveforms of the Lagrangian and the Hamiltonian. The Hamiltonian is fixed at 296.3 mJs with small parasitic variation, which is minimized via the proper choice of simulation options. The step ceiling is 1 ms.

The phase trajectories in Figure 12 and Figure 13 represent the contours of the generalized energy (47) and Hamiltonian (48). They can be discussed the same way as the trajectories in Figure 8 and Figure 9.

6. Discussion

Hamilton’s variational principle, i.e., that the trajectory is the extremal of Lagrange’s function, holds for the Σ-circuits [16]. This work offers new knowledge, namely that a Hamiltonian can be constructed for Σ-circuits and that this Hamiltonian preserves its fixed value during motion. Its physical dimension is given by the number Σ of the diagonal, and its unit is [VAs^—Σ]. The Hamiltonian of a circuit consisting of generally nonlinear inductors and capacitors, therefore, represents the energy because both the inductor and capacitor are located on the diagonal with the number Σ = −1. The physical character of the Hamiltonian does not change even if other types of elements from this diagonal appear in the circuit.

The work in [9] was the beginning of this type of research. This study introduced the state functions of memristors, memcapacitors, and meminductors and constructed a Lagrangian from those functions. The research in [11,12,16], which appeared later, introduced the conditions for validating Hamilton’s variational principle, building up Hamilton’s formalism, and generalizing all these fundamental pieces of knowledge for general higher-order elements. It was demonstrated in [16] that the higher-order Lagrangian is “native” for general HOEs. Our work offers another insight: For this Lagrangian we construct the corresponding Hamiltonian, which is “native” for circuits with HOEs.

In order to derive benefits from Hamilton’s formalism in circuits with HOEs, we should use Ostrogradsky’s formalism because only the higher-order Hamiltonian is native for such circuits. One of its advantages is the ability to construct comfortable analyses of the preserved quantities when finding periodical steady states or the conditions of the occurrence of self-oscillation, as demonstrated by the example of the PU oscillator. The synthesis of this circuit via nonlinear HOEs leads to a nonlinear oscillator, which has not been hitherto described in the literature. The Hamiltonian of the PU oscillator consisting of the FDNR, R, and FDNC elements, according to Figure 4a, represents the power because these elements are located on the diagonal with Σ = 0. The Hamiltonian of the PU oscillator consisting of MR, linear FDNC, and linear FDPC elements, according to Figure 4b, represents the action (integral of energy) because these elements are located on the diagonal with Σ = −2. This quantity is also preserved in the system in spite of the fact that the system is dissipative. Quod nota, the paradox is only illusory because two of the three elements of these specific circuits are, in principle, active elements, namely the FDNR and FDNC or the linear FDNC and linear FDPC elements.

The Hamiltonian of the circuit, containing only the HOEs from one diagonal, can be constructed via Equation (24). Equation (24) represents the known structure of a higher-order Hamiltonian, in which the generalized coordinates, generalized momenta, and higher-order Lagrangian appear. The Lagrangian is set up via the constitutive relations of the individual elements. The individual momenta are arranged according to (24) as sums of the generalized voltages of the corresponding orders across selected elements of the circuit. The Lagrangian is given as a sum of the state functions of the individual elements provided with the appropriate sign. The state function of the element corresponds to the area below its constitutive relation.

It has been newly found that the state function of the well-known HP memristor with the Joglekar window function for p = 1 is a composition of the dilogarithm and the exponential function.

The modeling diagram of the one-dimensional Σ-circuit from Figure 3 can be generalized in a classical way to a multi-dimensional system (e.g., a circuit with several loops), which can be used for implementation in the corresponding simulation program. For simulation, it is desirable to select a method for numerical integration that is suitable for the analysis of the sets of Hamilton’s differential equations [27]. The waveform of the Hamiltonian must exhibit negligible deflections from the constant value. This criterion provides a good feedback when seeking for an optimum configuration of the simulation task and parameters, including the time step.

The theoretical apparatus built in this work can also be used for circuits conisting only of two types of HOEs immediately neighboring each other on the common Σ-diagonal. However, the Hamiltonian of such a circuit is of a classical not a higher-order type. A typical example is a network of mutually interconnected memcapacitors and meminductors.

For researchers dealing with concrete memristive devices, for example the TiO_x, TaO_x, SiO_x, HfO_x and other types, it will be useful to determine whether the presented Hamilton’s formalism for higher-order elements can also be applied to these devices.

It follows from the essence of the given formalism that this formalism can be used only for circuits that are made up of two-terminal HOEs, whereas all the HOEs from the circuit must be located on the common Σ-diagonal of Chua’s table. The question is whether the models of these memristive systems can be built from the above HOEs and whether all the other elements from the circuit are also located on the common diagonal.

Historically, the first model of the HP memristor with a simple window function [7,28] for modeling the nonlinear dopant drift can be classified as an ideal generic memristor. Since this memristor is equivalent to an ideal memristor, it can be considered an HOE of the (−1,−1) type, which complies with Hamilton’s formalism. However, the corresponding application circuit must contain, in addition to this memristor, only the other appropriate HOEs from the Σ = −2 diagonal—for example, (−2, 0), (0, −2) and other elements. If some of these elements are linear, then they can be moved in Chua’s table along the corresponding Δ-diagonals [31]. Then, the set of these admissible linear elements will grow, as is illustrated in Section 5.

However, the complex physical models of the above mentioned memristive devices, such as the Tunneling barrier model, TEAM, and others [32], are classified as extended memristors, which cannot be generally built from two-terminal HOEs. Thus, attempts to utilize the presented Hamilton’s formalism for such systems interfere with the principal limits. When looking for a Hamilton’s formalism suitable for these elements, it will be necessary to use other methods, such as those mentioned, for example, in [21,22].