Para-Hamiltonian form for General Autonomous ODE Systems: Introductory Results

Abstract

We propose a new tool to deal with autonomous ODE systems for which the solution to the Hamiltonian inverse problem is not available in the usual, classical sense. Our approach allows a class of formally conserved quantities to be constructed for dynamical systems showing dissipative behavior and other, more general, phenomena. The only ingredients of this new framework are Hamiltonian geometric mechanics (to sustain certain desirable properties) and the direct reformulation of the notion of the derivative along the phase curve. This seemingly odd and inconsistent marriage of apparently remote ideas leads to the existence of the generator of motion for every autonomous ODE system. Having constructed the generator, we obtained the Lie invariance of the symplectic form $\omega$ for free. Various examples are presented, ranging from mathematics, classical mechanics, and thermodynamics, to chemical kinetics and population dynamics in biology. Applications of these ideas to geometric integration techniques of numerical analysis are suggested.

1. Introduction

Physics clearly distinguishes systems with conservation laws obeyed, as this simplifies the treatment of physical input given by initial or boundary conditions and a formal description in terms of admissible mathematical tools. Hamiltonian systems, which are model examples of how conserved quantities fit into physics, are very convenient since their nice formal features and associated structure provided by symplectic geometry additionally reveal several hidden links between objects engaged in the theory [1]. Nonetheless, often, the analytic methods are not sufficient, as we need to ponder the possibility of the numerical solution of the problem under consideration or even some trickier purely formal methods, as asymptotic expansion in the vicinity of singularities [2].

This article proposes some form of these “trickier methods” that could be applied to numerical algorithms, as they offer a weak form of the conservation laws (we mean here a conserved charge without the corresponding classical conservation law.

Numerical integration algorithms range from general-purpose classic Runge–Kutta and multi-step schemes [3] to more specific, geometric methods [4]. The latter admit of course better qualitative behavior (e.g., long time stability [5]), while the former usually are of higher order with less-flexible adaptation possibilities when it comes to the change of the simulation circumstances (e.g., the scheme fails to maintain high accuracy when the simulation lasts too long). Hence, geometric methods show some superiority, when we go to the extremes of bending the simulation parameters.

We found the main motivation of the undertaken research in the geometric numerical treatment of ODEs known as Geometric Numerical Integration (GNI) (see, e.g., [4,6,7,8,9]), although in this paper, we confine ourselves mainly to the construction of the proper continuous counterpart of the discrete framework needed for the treatment of general autonomous, non-conservative systems.

As is well known, geometric methods demand some kind of in-built qualitative structure of the considered problem to be preserved by the scheme. A good example of such a structure, and maybe the simplest one, is a single conserved quantity. This case can be reduced to the usage of a certain discrete gradient algorithm; see, e.g., [10,11,12,13]. The last reference here is an excellent example of accessing the second invariant quantity through discrete gradients. Such a numerical treatment was indeed applied to systems with first integrals and Lyapunov functions, but for arbitrary systems, it generally ceased to function properly because of a lack of a structural property guiding the evolution of the system. Although several frameworks of interest have been proposed for the qualitative description of such problems (see, e.g., [14,15,16,17,18]), none of these provides a simple and immediate way to gain access to the conserved quantity, helping directly in the performance of a specially crafted numerical algorithm.

In this paper, we show that Hamiltonian systems are indeed ubiquitous when we consider some consequences of the theorem on the existence and uniqueness of solutions to ordinary differential equations (ODEs) [19,20,21] and seek for the structure described above even if it is apparently absent. This approach leads to an abundance of so-called effective integrals of motion.

To begin with, let us consider the simple Initial-Value Problem (IVP) for an autonomous ODE:

$$\dot{\mathit{x}}=\mathit{f}\left(\mathit{x}\right),\mathit{x}\left(t_0\right)={\mathit{x}}_0,$$

(1)

given on some open domain $D\subset {\mathbb{R}}^d$ with $\mathit{f}:D\to {\mathbb{R}}^d$; obviously, an overset dot stands for a shortcoming denoting time derivation.

Given $\mathit{f}\left(\mathit{x}\right)\in {\mathcal{C}}^0\left(D\right)$, the existence and continuity of the solutions are guaranteed. The assumption of $\mathit{f}\left(\mathit{x}\right)\in {\mathcal{C}}^r\left(D\right),r⩾1$ ensures the uniqueness of the solution and its respective differentiability properties [1,19,21].

McLachlan et al. [12,22] showed that, in a neighborhood of a non-degenerate fixed point of a dynamical system, the existence of the first integral of (1) is equivalent to the existence and boundedness of skew-symmetric matrix $B\left(\mathit{x}\right)$ such that:

$$\mathit{f}\left(\mathit{x}\right)=B\left(\mathit{x}\right)\nabla H\left(\mathit{x}\right),$$

(2)

where $H\left(\mathit{x}\right)$ denotes the mentioned integral. Note that B is not determined uniquely, since we can add any solution of the homogeneous equation $0=B\nabla H$ to the particular solution:

$$S^{\left(P\right)}=\frac{1}{{|\nabla H|}^2}\mathit{f}\wedge \nabla H$$

(3)

The exterior product could be defined as:

$$v\wedge w:=\sum _{i,j}(v_i{w}_j-v_j{w}_i),i,j=1,2,\dots ,d,$$

(4)

being explicitly anti-symmetric (note that in ${\mathbb{R}}^3$, the structure provided by such a product is equivalent to the one given by the usual cross-product). The Formula (3) explicitly demands $|\nabla H|\ne 0$. However, if we are already given the problem in the form of (2), the assumption of $\nabla H\ne 0$ is redundant.

This article is meant to serve to extend the classical Hamiltonian approach to non-conservative systems, therefore unifying the perspective of gradient discretization for all autonomous ODEs and their study, as this has its origin in [22]. The approach dedicated to even more general systems (non-autonomous, optimal control, stochastic) is still under development.

In Section 2, we provide an overview of the subject of geometric mechanics (for the notation and basic terminology, check [1,23]), Section 3 serves in the same way to refresh the basic information on the biological population and chemical reaction dynamics. Section 4 introduces the main ideas of the paper, emphasizing direct overuse of the existence and uniqueness theorem for ODE systems. Section 5 presents versatile examples of the application of the introduced formalism, and Section 6 gives concluding remarks.

2. Hamiltonian Mechanics

If $B\left(\mathit{x}\right)$ in (2) obeys the Jacobi identity, we can term the system Poisson or Hamiltonian (non-canonical). Since Poisson structure matrix B is always of an even rank [24], odd-dimensional systems are necessarily degenerate (some even-dimensional ones also are; this implies the existence of Casimir functions) [25,26].

On the chosen symplectic leaf with prescribed Clebsch coordinates, treated as symplectic manifold $(T^{*}M,\omega )$, we have:

$$\omega =d{x}^i\wedge d{p}_i$$

(5)

(the summation convention assumed throughout the paper) the existence of which is equivalent to the non-degenerate canonical Poisson bracket given by the Poisson bi-vector:

$$\pi =-{\omega }^{-1}:\stackrel{2}{\bigwedge }{T}^{*}M\to \mathbb{R}.$$

(6)

where:

$$\pi (df,dg)=\{f,g\},$$

(7)

and the $f,g\in \mathcal{F}\left(T^{*}M\right)$ Lie algebra of functions defined on the phase space.

This also gives rise to the first-order differential operator:

$$X_H(·)=\pi (·,dH)$$

(8)

known as a Hamiltonian vector field.

Now, we call it the the first integral generator of motion if:

$$\dot{\mathit{x}}=\{\mathit{x},H\}=X_H\left(\mathit{x}\right)$$

(9)

understood componentwise. In other words:

$${\dot{x}}^i=\frac{\partial H}{\partial p_i},{\dot{p}}_i=-\frac{\partial H}{\partial x^i}.$$

(10)

or simply:

$$X_H\rfloor \omega =dH,$$

(11)

denoting by ⌋ the substitution of a vector field into the form $\omega$ (contraction).

Note that the flow of the Hamiltonian vector field preserves the canonical symplectic form on the phase space, which is clearly given by the proper Lie–Ślebodziński derivative:

$${\mathcal{L}}_{X_H}\left(\omega \right)=X_H\rfloor d\omega +d(X_H\rfloor \omega )=d\left(dH\right)=0,$$

(12)

by the closedness of the symplectic form and the nilpotency of exterior derivative d.

The basic hydrodynamic interpretation of the canonical formalism is also of some value. Let us consider a phase fluid of many systems with various initial conditions, moving on the phase space [27]. The velocity field of the fluid is clearly given by the Hamiltonian vector field. Note that since $H\in {\mathcal{C}}^2\left(D\right)$, we have $\mathrm{div}\mathit{v}=0$.

It surely obeys the continuity law:

$$\frac{\partial \rho }{\partial t}+\mathrm{div}\left(\rho \mathit{v}\right)=0,$$

(13)

or, in other words:

$$\frac{d\rho }{dt}+\rho (\nabla ·\mathit{v})=0.$$

(14)

Now, because the phase fluid is incompressible:

$$\rho =\mathrm{const}$$

(15)

on a target energy level set. Now, let us ponder $\rho =C$, and then, we vary the constant as $C=C(x,p)$ to obtain the condition:

$$\nabla C·\mathit{v}=0\to C=cH,$$

(16)

where c is a rightful constant.

The fluid should also undergo Euler’s equation (the fluid’s particle velocity denoted by $\mathit{v}$; here, we considered just the fluid of an ensemble of free harmonic oscillators with the equations of motion $\dot{x}=p,\dot{p}=-x$):

$$\frac{d\mathit{v}}{dt}=-\frac{1}{\rho }\nabla P$$

(17)

which in an elementary manner leads to Bernoulli’s law, if we take the inner product of both sides with $\mathit{v}$ and perform trivial integration:

$$\frac{1}{2}c(x^2+p^2)+P\left(\mathit{x}\right)=\mathrm{const},$$

(18)

from which we can evaluate the pressure function and constant appearing in the formula, as the pressure needs to be non-negative.

3. $\mathcal{M}$-Systems of Ecology and Chemical Kinetics

The Hamiltonian (in general Poisson) structure can be met also in ecology, where evolving populations of different species share resources in a common domain of living. Various environmental factors also can be modeled through some additional terms in differential equations governing the evolution (optimal control problems, etc. [28]).

We point out that a similar approach can be adopted for chemical kinetics problems, where the role of populations is played by the concentrations of various different chemical substances. The resource function here, if it exists, reflects the amount of reacting substances. In both domains, the systems share the properties of non-negativity, realizability, reducibility, and semi-stability [6,29,30].

The mentioned processes may involve a great number of variables, if they are to be described exactly. For our purposes, we stick to ODE models, as they are fair enough to describe phenomena with satisfying accuracy, yet simple enough not to complicate things unnecessarily, so we used the following assumptions:

1.

Reagents (species) are well mixed (distributed homogeneously), otherwise the problem would be inhomogeneous in space, hence yielding Partial Differential Equations (PDEs) (e.g., the reaction–diffusion problem; see, e.g., [31]) instead of ODEs;

2.

The concentrations of the substrates (species) are high enough to prevent stochastic behavior during reacting (coexisting; the abandonment of this assumption would lead to Wiener processes). Otherwise, we would obtain Stochastic Differential Equations (SDEs) (e.g., the Kubo oscillator) instead of ODEs.

Generally we could consider a few types of interaction between species: parasitic invasion, competition for resources, etc. For example, the competition of U and V given by the system:

$$\begin{array}{c}\dot{u}=f(u,v),\hfill \\ \dot{v}=g(u,v)\hfill \end{array}$$

(19)

would undergo conditions $f_u>0,g_v>0,f_v<0,g_u<0$ to reflect the fact that the consumption of a nutrient/the depletion of a resource by one species would prevent the other from doing the same, as well as describing the competition between the members of the same population.

Typically considered in this context is the Lotka–Volterra system of the form:

$$\begin{array}{c}\dot{u}=-\alpha u+\beta uv,\hfill \\ \dot{v}=-\gamma uv+\delta v,\hfill \end{array}$$

(20)

with a properly adjusted set of constants (note that the above equations fulfill the previously mentioned restrictions, if $u,v$ are sufficiently bounded) and u being the population of the predator species and the v-population concentration of the prey.

A specific kind of interaction appears in systems of the Rosenzweig–MacArthur [32,33] type, involving predator and prey coexistence. In (19), we would obtain:

$$\begin{array}{c}{\displaystyle f(u,v)=ru\left(1-\frac{u}{K}\right)-vh\left(u\right),}\hfill \\ {\displaystyle g(u,v)=v(-\beta +\alpha h(u\left)\right),}\hfill \end{array}$$

(21)

where $h\left(u\right)$ denotes the number of prey caught by the predator per unit of time and $\alpha ,\beta ,r,K$ are the constant environmental parameters of the system.

The described model may be decomposed into symmetric and anti-symmetric part-evolutions, although it does not give up the Poissonian description at all. The ambitious objective of this paper was to consequently remedy the situation, putting (21) into not only the Poissonian, but the altogether canonical form. This type of impossibility of the direct integration of the equations of motion to solve the inverse Hamiltonian problem renders the main motivation for all efforts considered here.

Functions as $h\left(u\right)$ are known as the functional response of the predator to a variation in the prey population. The basic types of these were classified by C.S. Holling into three categories: I. linear, II. hyperbolic (saturation), and III. so-called $\theta$-sigmoid [33].

Various parameters in the equations of the population dynamics can be made into variable ones, often leading to their re-appearance as new dynamical (dependent) variables. To study these and more of not only single-population models such as, e.g., a chemostat system with continuous/batch cultivated bacteria, the Monod equation, the grazers–vegetation cycles, the Ivlev or Ayala–Gilpin–Ehrenfeld model of population growth, epidemic and endemic (SIR/SIS) models, or even more complicated problems associated with the optimal control of invasive species and the harvesting of populations, the reader is referred to the literature on the subject [28,34,35].

Of course, some of these systems admit Hamiltonian/Poisson representation, where we treat the resource function $\mathcal{M}(\mathit{x},t)$ as the basic object (hence the name $\mathcal{M}$-systems). It serves as a generator for the equations of motion:

$$\dot{\mathit{x}}=B(\mathit{x},t)\nabla \mathcal{M}.$$

(22)

For instance, in the case of the Lotka-Volterra system (20) we have:

$$\mathit{B}=\left(\begin{array}{cc}0& uv\\ -uv& 0\end{array}\right)$$

(23)

with the generator:

$$\mathcal{M}(u,v)=-\alpha lnv-\delta lnu+\beta v+\gamma u$$

(24)

yielding the fine Poisson system or the Hamiltonian system of a non-canonical form. Sometimes, the term “M-systems” is used to describe molecular systems in biology; however, we stick to the meaning proposed in [36] and accepted by some other authors (e.g., [37]).

Remark 1.

Similarly, we can formulate the Hamiltonian-like $\mathcal{M}$-system substituting J in place of B and properly transforming the variables. However, a big difference occurs when we ponder both Poisson and Hamilton formulations of the same problem: the concentration variables for different species should be of non-negative values, but in the Hamiltonian case are not (incompressible fluid on ${\mathbb{R}}^m$). Because of that, we can interpret the compressibility of the phase fluid in the Poisson case as partially arising from constraining the phase space to be ${\mathbb{R}}_{+}^m$.

The set of chemical reactions is termed the reaction network. Having:

$$A_i+E_i\to P_i+B_i,i=1,\dots ,n,$$

(25)

we call species on the left the reactants and species on the right products of the reaction. When all stoichiometric coefficients are equal to one, we call such a reaction an elementary reaction. Note that every particular reaction can be written as an elementary reaction when we substitute $A_{1}X=X+\dots +X$ ($A_1$ times). As an example, we may consider the Robertson reaction network, mentioned in Section 5, and there should be no confusion in retaining the system in the potential, Poissonian form.

Chemical reactions formulated as a population dynamics problem use the mass action law [29,35]: at constant temperature, for any elementary reaction, its rate is proportional to the concentrations of the reactants.

The matrix formulation is also accessible to problems concerning the mass action law; see, e.g., [6,29]. This simple rule underlies the differential description of such reactions as given in the Michaelis–Menten model, Hill enzymatic equation, or Robertson network.

One of the fundamental features of mass action kinetics is that it produces differential equations with polynomial non-linearities. This also means that when we encounter such a set of equations, we may find the reaction network obeying these. Such a process is referred to as the realizability of the mass action kinetics [29]. An example of this procedure may be the so-called Lotka–Volterra reactions, retrieved from (76); see [6,38].

4. Para-Hamiltonian Description of Non-Conservative Systems

Basically speaking, we pondered a classical system with energy $E=\frac{1}{2}{p}^2+V\left(x\right)$, where the number of dimensions is irrelevant (for now). However, in terms of Newtonian forces, in general, they are essentially non-self-adjoint (simply meaning these are not integrable [39]):

$$\dot{x}=p,\dot{p}=-\mathrm{grad}V\left(x\right)-D(x,p),$$

(26)

where D stands for dissipative forces.

The theory of ODEs confirms that if we are able to find a unique solution to the problem, then the phase trajectories do not intersect, in order to keep the vector field generating the ODE well defined. This means, provided sufficiently differentiable V and D, that any given instant t in time connects by the 1:1 correspondence to some $(x,p)$, and vice versa. The following exposition will lean heavily upon this fact and exploit the mentioned correspondences in terms of the particular differentiation procedure along the trajectory. As an effect, we construct a new class of effectively conserved quantities. Hence, we assumed that $V\in {\mathcal{C}}^{\infty },D\in {\mathcal{C}}^{\infty }$ in all their arguments in a given domain, in a classical sense.

Definition 1.

Given a Pfaff form $\theta =\kappa (x,p)dx+\nu (x,p)dp$ (where κ and ν are some functions) and a system of equations whose trajectory (parameterized by t) maps initial data $(x_0,p_0)$ to $(x,p)$, we define:

$$\Psi (x,p)=\Psi (x_0,p_0)+{\int }_{\gamma (t_0,t)}\theta \equiv \Psi (x_0,p_0)+{\int }_{(x_0,p_0)}^{(x,p)}\left(\kappa (x^{\prime },p^{\prime })d{x}^{\prime }+\nu (x^{\prime },p^{\prime })d{p}^{\prime }\right),$$

(27)

where the path of the integration is along this trajectory, denoted by $\gamma (t_0,t)$. Then, we denote $đ\Psi :=\theta$, and:

$$\frac{\partial \Psi }{\partial x\left(t\right)}:=\kappa (x\left(t\right),p\left(t\right)),\frac{\partial \Psi }{\partial p\left(t\right)}:=\nu (x\left(t\right),p\left(t\right)).$$

(28)

Note that the above expressions are not standard partial derivatives, except in the case when the Pfaff form θ is an exact differential. In the following, these operations are referred to as derivatives along the trajectory.

We start with the equations of motion (26) as given. Then, (assuming, for simplicity, the one-dimensional case) we seek for a “generator of motion”, however, not in a standard sense, but in the sense of Definition 1 (this is what me mean by the “para-Hamiltonian” description).

$$\begin{array}{c}{\displaystyle \dot{x}=p\stackrel{?}{=}\frac{\partial K}{\partial p\left(t\right)}\to K=\frac{1}{2}{p}^2+c_1\left(x\right),}\hfill \\ {\displaystyle \dot{p}=-V^{\prime }\left(x\right)-D(x,p)\stackrel{?}{=}-\frac{\partial K}{\partial x\left(t\right)}\to K=V\left(x\right)+{\int }_{\gamma (t_0,t)}D(x,p)dx+c_2\left(p\right),}\hfill \end{array}$$

(29)

where $c_1\left(x\right)$ and $c_2\left(p\right)$ can be derived by the comparison of both equations. Hence:

$$K=\frac{1}{2}{p}^2+V\left(x\right)+{\int }_{\gamma (t_0,t)}D(x,p)dx$$

(30)

In other words, we define in the phase space $(x,p)$ one-form $đK$, which is neither closed nor exact:

$$đK=\frac{\partial K}{\partial x\left(t\right)}dx+\frac{\partial K}{\partial p\left(t\right)}dp=(V^{\prime }\left(x\right)+D(x,p))dx+pdp.$$

(31)

There is ongoing energy exchange with the environment (not necessarily a loss) due to non-potential function $D(x,p)$. We introduce a weird object w called the “reservoir”, dependent on x through the upper limit of integration. Its x-derivative along the trajectory is $D(x,p)$, and theoretically, it is a redundant variable; hence:

$$w={\int }_{x_0}^{x\left(t\right)}D(x^{\prime }\left(t^{\prime }\right),p^{\prime }\left(t^{\prime }\left(x^{\prime }\right)\right))d{x}^{\prime }\left(t^{\prime }\right),$$

(32)

where all needed inter-dependencies are guaranteed by the existence and uniqueness theorem for ODEs [1]. In the above formula, we perform the Riemann–Stieltjes integral to guarantee there exists some form of 1:1 correspondence between dynamical variables.

In practice, adjoining the reservoir to the system means adjoining non-potential (dissipative) forces as working for the positive account of the energy of the system (meaning they are treated as part of the system). By means of physical analysis, we incorporate the “power continuity law”:

$$\frac{dE}{dt}+pD(x,p)=0,$$

(33)

also giving rise to one-form $đK$, which is perceived as a fundamental quantity of our approach (consequently, we vary only the dependent variables).

Meanwhile, it is worth noting that among the many obtained result aimed at grasping dissipative behavior correctly, the Nosé–Hoover system [13] and its generalizations [40] deserve some amount of attention. While cast in a Hamiltonian-like form (a not-so-canonical Poisson tensor), the equations of motion link the dynamics of the system to the state of the thermostat (through temperature T), which gives effective access to the thermodynamics of the ensemble of such systems (although not the statistical dynamics of such, the divergence of the vector field is in general non-zero). The importance of this remark rests on the similarities of this description with the one proposed here.

Theorem 1.

The quantity:

$$K=E+w$$

(34)

is conserved.

Proof. 

We may easily check:

$$\dot{K}=\dot{E}+\dot{w}=-D(x,p)p+D(x,p)p=0;$$

(35)

however, it is not a well-defined function of x and p, since it depends on the path taken by the system. □

The above result is important not because of its complicated nature, but because of its simplicity. We accessed a novel type of conserved quantity that is rather of no use in pure theoretical considerations. As outlined in the Introduction, it is of huge practical/computational value.

Let us observe that K possesses an extremely trivial physical interpretation: it is just the initial energy (provided that $w\left(t_0\right)=0$). Since, for general $(x_0,p_0),w\left(t_0\right)=\mathrm{const}$, we have $K=E(x_0,p_0)+w\left(t_0\right)$, therefore, K is a smooth function of initial values provided that E is smooth.

We call the differential form not being the differential of any function a (pure) Pfaffian form [41]. An example of such a quantity is:

$$đw=D(x,p)dx.$$

(36)

In the future, we will denote by $w_x$ the expression standing next to $dx$ in $đw$.

Therefore, we can define:

$$đK=đE+đw,$$

(37)

as a Pfaffian differential form.

In the above definition, we paid attention to the fact that the Hamiltonian is no longer preserved. Its decrement is exactly the increment of w with the opposite sign. Since the change of the Hamiltonian now obviously depends on the path taken by the system on the phase space, we cannot even claim that the Hamiltonian is still a potential-type function or a properly defined function at all, although, as a shortcoming, we use the term “potential part of the generator of motion” with respect to the Hamiltonian (so, strictly speaking, $dE$ becomes the Pfaffian form).

Since $\dot{E}=-D(x,p)p$ (see (33)), the addition of a reservoir with the exactly opposite time derivative yields the conserved quantity. Consequently, we tried to describe the problem in two different sets of canonical coordinates in several dimensions: $(x^i,p_i),(u^i,w_i),i=1,\dots ,n$. Motion begins with $K_0=L_0=E(x_0,p_0)$; energy as a number at a given time is invariant with respect to the change of coordinates (“frozen” time transformations); since the initial value $K_0=L_0$ is also invariant, it induces the invariance of the integral $w={\int }_{x_0}^x{D}_i(x,p)d{x}^i$ while changing the description of the system so that:

$$\begin{array}{c}{\displaystyle đK=p_{i}d{p}_i+{\nabla }_{x^i}V\left(x\right)d{x}^i+D_i(x,p)d{x}^i,}\hfill \\ {\displaystyle đL=w_{i}d{w}_i+{\nabla }_{u^i}\tilde{V}\left(u\right)d{u}^i+F_i(u,w)d{u}^i,}\hfill \end{array}$$

(38)

where the summation over repeating indices is assumed.

Demanding the para-Hamiltonian framework to hold in both sets of symplectic coordinates:

$$\begin{array}{cc}{\displaystyle {\dot{x}}^i=\frac{\partial K}{\partial p_i\left(t\right)},}\hfill & {\displaystyle {\dot{u}}^i=\frac{\partial L}{\partial w_i\left(t\right)},}\hfill \\ {\displaystyle {\dot{p}}_i=-\frac{\partial K}{\partial x^i\left(t\right)},}\hfill & {\displaystyle {\dot{w}}_i=-\frac{\partial L}{\partial u^i\left(t\right)}}\hfill \end{array}$$

(39)

we obtain the condition for action differentials to differ by an exact differential (as usual):

$$p_{i}d{x}^i-w_{i}d{u}^i=d\varphi ,$$

(40)

where the time-dependent differentials cancel due to the common initial conditions. Traditionally performing proper Legendre transforms, we can express $\varphi$ as a function of the preferred two subsets of the old/new variables.

Here, the integral parts with the non-potential functions match due to the invariant character of energy; the covectorial character of momenta induces $w_j=p_i(x,w)\frac{\partial p_i}{\partial w_j}$; the left-over condition gives:

$$\left(D_i(x,p)-F_j(u,w))\frac{\partial u^j}{\partial x^i}\right)d{x}^i=0,$$

(41)

where we tacitly assumed $u\left(x_0\right)=u_0$. The demand for the integral to vanish is sufficient for the integrand to vanish as well. The vanishing of the integral demand is sufficient for the integrand to also vanish. However, the dependence of D on p (or w) forbids its form as a simple gradient function. Hence, D is a covector, but it is not a gradient of any scalar function.

The potential part is invariant due to the symplectic structure; hence, the dissipative form integrated from $x_0$ to x or from $u\left(x_0\right)$ to $u\left(x\right)$ has to be invariant in order to ensure the independence of the energy dissipated during the motion. We call such invariance secondary (or induced) invariance.

From earlier considerations, we may generalize (8) to:

$$X_K=\pi (·,đK)$$

(42)

giving rise to the Poisson bracket:

$$\{f,K\}=X_K\left(f\right),$$

(43)

for any function given on the phase space. We can conceive of K as the formal generator of non-potential motion: its Poisson bracket with canonical coordinates gives proper equations of motion:

$$\begin{array}{c}\dot{x}=\{x,K\}=\pi (dx,đK)=p,\hfill \\ \dot{p}=\{p,K\}=\pi (dp,đK)=-V^{\prime }\left(x\right)-D(x,p)\hfill \end{array}$$

(44)

moreover:

$$\dot{w}=\{w,K\}=\pi (đw,đK)=D(x,p)p,$$

(45)

which is indeed the case.

We can regard the condition:

$$đK\left(\mathit{v}\right)=0$$

(46)

for any $\mathit{v}\in T\left(T^{*}M\right)$ as a formal substitute for the law of the conservation of energy, determining the trajectory uniquely.

In agreement with the thermodynamic description, we may think of the energy function as an analogue of the internal energy U. Thus, we can yield two different interpretations:

.1

Adjoining external forces to the language in which we describe the system, we obtain a closed system; hence, $dU=-PdV=-đw=D(x,p)dx$, and hence, the non-potential force exerts a kind of “pressure” on the ensemble of systems in the phase space (see below);

2.

If $đdK=dU=0$, then we view $TdS$ and $PdV$ as equal, so reaching for the phase pressure obtained below for the phase fluid and considering phase volume $dxdp$ as $dV$, we can pick the monotonically growing component of the expression $PdV$ and treat it as an infinitesimal increment of the entropy analogue of the system, the remaining factor being the “temperature”.

Summing up all the differentiability conditions, we seek $E(x,p)\in {\mathcal{C}}^n,w\left(x\right)\in {\left(C\right)}^1$ such that $K=E+w\in {\mathcal{C}}^1$, but $\nabla K\in {\mathcal{C}}^m$, $m,n$ arbitrary.

With this setting in mind, we can derive new formulas for the vector field algebra. Considering K the formal generator of motion and f a well-differentiable function, we obtain:

$${[X_f,X_K]}_L=-X_{\{f,K\}}+\mathrm{div}\mathit{v}{X}_f,$$

(47)

seeing that compressible terms ($\mathit{v}$ is of course a tangent vector to the phase flow, namely $X_K$) that are responsible for the discontinuities are causing an anomaly in the vector field algebra to occur.

It is possible to find a similar formula for the pair of reservoir-containing $K,L$, say. However, it is not very useful in the context of the Hamiltonian description of mechanical systems since there is only one of these needed to govern the dynamics. The situation dramatically changes in Nambu, or generalized Nambu, mechanics (e.g., [42]), where the dynamics is given in terms of a few such vector fields.

We stick to this working approach, especially since it guarantees:

$$\omega =dx\wedge dp$$

(48)

as a symplectic form. During the former considerations, it was preserved by the flow of Hamiltonian vector field $X_H$, and now, it satisfies:

$${\mathcal{L}}_{X_K}\left(\omega \right)=X_K\rfloor d\omega +d(X_K\rfloor \omega )=d(đE+đw)=d(-đw+đw)=0.$$

(49)

as it provided by the equations of motion: $\dot{p}=-V^{\prime }\left(x\right)-D(x,p)$, hence $dp=-V{;}^{\prime }\left(x\right)dt-D(x,p)dt$, so multiplying both sides by p, we obtain $đE=pdp+V^{\prime }\left(x\right)dx=-D(x,p)dx$.

$X_K$ is defined unambiguously throughout the phase space as a section over $T\left(T^{*}M\right)$, hence providing that phase trajectories do not cross each other.

When it comes to the value of K, it is determined by providing initial conditions $({\mathit{q}}_0,{\mathit{p}}_0)$. Hence, $K\in \mathcal{F}\left(T^{*}M\left(t_0\right)\right)$, so its constant value is uniquely determined by the state of the system at the initial moment. As a function of the flow, K may be given as:

$$K={\int }_{t_0}^t\mathit{v}\rfloor \omega$$

(50)

making its dependence on the initial conditions much less manifest.

A little bit more sophisticated is the simultaneous use of two reservoirs for the system:

$$\begin{array}{c}\dot{x}=p+F(x,p),\hfill \\ \dot{p}=-V^{\prime }\left(x\right)-D(x,p),\hfill \end{array}$$

(51)

yielding:

$$đK=(p+F(x,p))dp+(V^{\prime }\left(x\right)+D(x,p))dx$$

(52)

and, accordingly:

$$X_K=(p+F(x,p))\frac{\partial }{\partial x}-(V^{\prime }\left(x\right)+D(x,p))\frac{\partial }{\partial p}.$$

(53)

Since we are playing with ODE system, when we assume that $A\left(\mathit{x}\right)\nabla K=\mathit{f}\left(\mathit{x}\right)\in {\mathcal{C}}^r\left(D\right),r⩾1$, the theorem on the existence and uniqueness of solutions is in power [1,19,20]. Therefore, there is a 1:1 correspondence between every moment in time and the points in the phase space. Trajectories on the phase space of the system are obviously not crossing. Hence, we can perform in an unambiguous sense any integral of a function of variables of the system with respect to some of these variables (or time) as a Riemann–Stieltjes integral.

Therefore, we can write:

$$K=\frac{1}{2}{p}^2+V\left(x\right)+w+z,$$

(54)

where the reservoirs are defined by:

$$\begin{array}{c}{\displaystyle w={\int }_{t_0}^{t}D(x\left(t^{\prime }\right),p\left(t^{\prime }\right))dx\left(t^{\prime }\right),}\hfill \\ {\displaystyle z={\int }_{t_0}^{t}F(x\left(t^{\prime }\right),p\left(t^{\prime }\right))dp\left(t^{\prime }\right),}\hfill \end{array}$$

(55)

so that $\dot{K}=0$.

Now, in order to make the current discussion as similar to the conservative case as possible, we focus for a moment on the hydrodynamical analogy, starting from the continuity equation:

$$\frac{d\rho }{dt}+\rho \nabla ·\mathit{v}=0,$$

(56)

hence:

$$\rho =C{e}^{{\int }_{t_0}^t{D}_p(q,p)d{t}^{\prime }}$$

(57)

where C is a constant and the lower index denotes the derivative with respect to an argument.

The integral in the exponent does not cause any trouble, since all the fluid’s particles obey the equations of motion; hence, we may again apply the Riemann–Stieltjes integral.

Note that we can consider C as depending on the canonical variables, where from the continuity equation, we obtain the constraint $\mathit{v}·\nabla C=0$, so C may depend on the K value, so on the phase trajectory of interest. Here, as always, $\mathit{v}=X_K$.

Additionally, we have Bernoulli’s law (from Euler’s equation):

$$\frac{1}{2}\rho {\mathit{v}}^2+P=\mathrm{const}.$$

(58)

Taking an example of linearly damped harmonic oscillators with the equation of motion:

$$\begin{array}{c}\dot{x}=p,\hfill \\ \dot{p}=-x-bp.\hfill \end{array}$$

(59)

Since all fluid particles obey these equations of motion, the continuity equation yields:

$$\rho =cK{e}^{bt},$$

(60)

where c is a constant.

Writing Bernoulli’s law:

$$\frac{1}{2}cK{e}^{bt}(x^2+p^2+2bpx+b^2{p}^2)+P=\mathrm{const},$$

(61)

and remembering that $E\sim e^{-bt}$, we see:

$$P=P_0-cK{e}^{bt}(bxp+\frac{1}{2}{b}^2{p}^2).$$

(62)

Fortunately, we know solutions to the damped oscillator, being $x\sim e^{-\frac{b}{2}t}(cos(\omega t+\delta ))$; hence, $p\sim e^{-\frac{b}{2}t}(cos(\omega t+\delta )-\omega sin(\omega t+\delta ))$. Thus, we that see there is no danger of the variable part of pressure growing to infinity. Provided that the engaged constants obey:

$$P_0-cK{A}_0^2\left((b+\frac{b^2}{2}){cos}^2\delta -(\omega b+b^2)sin\delta cos\delta +\frac{1}{2}{b}^2{\omega }^2{sin}^2\delta \right)⩾0,$$

(63)

where $A_0$ is the initial amplitude, the pressure is always positive. Notice that for different K (initial energy), this demand can somewhat change quantitatively.

5. Particular Non-Potential Systems

Here, we give a sequence of illustrative examples, treatable along the lines of the presented approach. Their objective is to show the applicability of the invented framework to low-dimensional systems and to directly show the presence and surprising preservation of the symplectic structure. As mentioned in the Introduction, the general statement with the proof will be published elsewhere.

Example 1. van der Pol oscillator

 

The van der Pol oscillator is a system that arises as some generalization of an RLC circuit (through the so-called Liénard-form equation of the VdP oscillator [19]), and its equations of motion are:

$$\begin{array}{c}\dot{x}=p,\hfill \\ \dot{p}=-x+\epsilon (1-x^2)p.\hfill \end{array}$$

(64)

The non-potential generator of motion is:

$$K=\frac{1}{2}(x^2+p^2)+w,$$

(65)

where the reservoir variable is given by:

$$w=-\epsilon {\int }_{t_0}^t(1-x^2\left(t^{\prime }\right))p\left(t^{\prime }\right)dx\left(t^{\prime }\right),$$

(66)

turning K into an effectively conserved quantity.

Remark 2.

Note that taking the $\partial K\left(t\right)$ derivative along the phase curve, we have:

$$\left(\begin{array}{c}\dot{x}\hfill \\ \dot{p}\hfill \end{array}\right)=\left(\begin{array}{c}\frac{\partial K}{\partial p\left(t\right)}\hfill \\ -\frac{\partial K}{\partial x\left(t\right)}\hfill \end{array}\right)=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)\nabla K$$

(67)

so that the system possesses a symplectic structure. Nevertheless, it may seem as it should not hold along the trajectory in the phase space, since $K\notin {\mathcal{C}}^2\left(T^{*}M\right)$ (exactly, we have $K\in {\mathcal{C}}^1\left(T^{*}M\right)$). However, due to the vanishing of $đK$ by the construction, we have (49), and the symplectic structure is preserved.

Remark 3.

Note that onward, to avoid cumbersome writing in integrals such as (66), we often write $y^{\prime }$ instead of $y\left(t^{\prime }\right)$, etc.

Example 2. Brusselator

 

The Brusselator is the dynamical system modeling the auto-catalytic reaction network [21]:

$$A\to X,2X+Y\to 3X,B+X\to Y+D,X\to E.$$

(68)

We claim that substrates $A,B$ are abundant in the environment, so we can denote their concentrations $a,b$ as being constant. We identify the X and Y species’ concentrations as dynamical variables x and y, respectively. From (68), with the use of the mass action principle, we obtain:

$$\begin{array}{c}\dot{x}=a+x^{2}y-bx-x,\hfill \\ \dot{y}=bx-x^{2}y.\hfill \end{array}$$

(69)

A quick look at these confirms there is no resource function $\mathcal{M}$ in the usual sense; the non-potential generator of motion becomes:

$$K=ay-\frac{1}{2}b{x}^2+w+z,$$

(70)

where the reservoir variables are given by:

$$\begin{array}{c}{\displaystyle w={\int }_{t_0}^{x\left(t\right)}{\left(x^{\prime }\right)}^2{y}^{\prime }d{x}^{\prime }}\hfill \\ {\displaystyle z={\int }_{t_0}^{y\left(t\right)}({\left(x^{\prime }\right)}^2{y}^{\prime }-b{x}^{\prime }-x^{\prime })d{y}^{\prime },}\hfill \end{array}$$

(71)

so that $\dot{K}=0$.

Remark 4.

Note that again, in the $x,y$ variables, the system is in canonical form $\dot{\mathit{x}}=J\nabla K\left(t\right)$ with:

$$J=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right),$$

(72)

which is preserved due to the courtesy of K being constant during the evolution.

Note that (69) has a fixed point at $(a,\frac{b}{a})$. This equilibrium is unstable when $b>1+a^2$; if $b<1+a^2$, it is stable. The case $b=1+a^2$ presents some doubts: in this situation, the origin appears to be the center (from the procedure of linearization); however, we know that if the dimension (here: the number of dependent variables) of the ODE system $n⩾2$, then the Hartman–Grobman theorem on linearization often fails at predicting the existence of the center [19,20,21].

Example 3. The Nosé–Hoover system [13,40]

We put the thermostatted system in the form:

$$\begin{array}{c}{\displaystyle \dot{q}=\frac{p}{m}=\frac{\partial K}{\partial p\left(t\right)},}\hfill \\ {\displaystyle \dot{p}=-{\Phi }_q-\frac{p{p}_{\eta }}{Q}=-\frac{\partial K}{\partial q\left(t\right)},}\hfill \\ {\displaystyle \dot{\eta }=\frac{p_{\eta }}{Q}=\frac{\partial K}{\partial p_{\eta }\left(t\right)},}\hfill \\ {\displaystyle {\dot{p}}_{\eta }=\frac{p^2}{m}-kT=-\frac{\partial K}{\partial \eta \left(t\right)},}\hfill \end{array}$$

(73)

which gives:

$$K=\frac{p^2}{2m}+\Phi \left(q\right)+\frac{p_{\eta }^2}{2Q}+kT\eta +{\int }_{q_0}^{q\left(t\right)}\frac{p^{\prime }{p}_{{\eta }^{\prime }}^{\prime }}{Q}d{q}^{\prime }-{\int }_{{\eta }_0}^{\eta \left(t\right)}\frac{p^{\prime 2}}{m}d{\eta }^{\prime },$$

(74)

from which, by differentiation, we can retain the equations of motion. Note that the transition to the time parameterization of the reservoirs (integral variables) shows simply $K=H$, the Hamiltonian presented, for example, by Ezra [13].

Remark 5.

The almost-Poisson structure provided in [40], given by:

$$B\left(\mathit{x}\right)=\left(\begin{array}{cccc}0& 0& 1& \frac{\tau p}{m}\\ 0& 0& 0& 1\\ -1& 0& 0& -p\\ -\frac{\tau p}{m}& -1& p& 0\end{array}\right)$$

(75)

is obviously of the non-canonical form. The time-dimensional parameter τ governing the evolution of the Nosé variable η must be set to zero, to obtain again the simple dynamics of (73).

Finally, we note that in our formulation of non-potential Hamiltonian mechanics, the system is canonical with two physical degrees of freedom. This property is obviously preserved along the evolution.

Example 4. Lotka–Volterra system

 

The Lotka–Volterra model describes the basic predator–prey interaction (with a linear response):

$$\begin{array}{c}\dot{u}=u-\alpha uv,\hfill \\ \dot{v}=-\beta v+uv,\hfill \end{array}$$

(76)

where u is the prey concentration in an environment, v is the predator concentration, α being the rate at which the consumption of the prey by a predator proceeds, and β is the death rate of a predator. Note that we chose the unit rate for the birth of the prey and the predator feeding on the prey.

We can write down these equations as a Poisson system:

$$\dot{\mathit{x}}=B\left(\mathit{x}\right)\nabla \mathcal{M},B\left(\mathit{x}\right)=\left(\begin{array}{cc}0\hfill & uv\hfill \\ -uv\hfill & 0\hfill \end{array}\right)$$

(77)

claiming that $\mathit{x}={(u,v)}^T,\nabla ={({\partial }_u,{\partial }_v)}^T$, and the resource function is:

$$\mathcal{M}=\beta lnu+lnv-u-\alpha v.$$

(78)

We should observe that such a formulated LV problem is given on the phase space ${\mathbb{R}}_{+}^2$ and, as a Poisson system, has a compressible phase fluid: $\mathrm{div}\dot{\mathit{x}}=1-\beta -\alpha v+u$; compare Remark 1.

We used the Poisson structure as even-dimensional and non-degenerate, so we can bring the system to its canonical form by transformation $u=e^q,v=e^p$, where the equations of motion become:

$$\begin{array}{c}\dot{q}=1-\alpha e^p,\hfill \\ \dot{p}=-\beta +e^q,\hfill \end{array}$$

(79)

with the incompressible phase fluid on the ${\mathbb{R}}^2$ symplectic phase space with the separable resource function:

$$\mathcal{M}=p-\alpha e^p+\beta q-e^q.$$

(80)

Note that the LV Poisson system would become of the canonical form also with:

$$đK=u(1-\alpha v)dv+v(\beta -u)du.$$

(81)

Moreover, we have:

$$\nabla K^{T}B\left(\mathit{x}\right)\nabla \mathcal{M}=0=\nabla {\mathcal{M}}^{T}J\nabla K;$$

(82)

hence, the integrals of motion commute in terms of each Poisson structure, but this only preserves the equilibria. What is more important is that we have the following.

Corollary 1.

The transition from the Poisson dynamics governed by the resource function $\mathcal{M}$ in (78) to the canonical form, the evolution of which is dictated by (81), preserves the Poisson bracket of a target function with the generator of motion (with its coupled matrix structure).

This remark is easily verifiable on a case-by-case basis, provided that the Poisson bracket of $\mathcal{M}$ with coordinates $u,v$ is preserved.

Knowing that the non-potential generator K provokes the anomalies of the vector fields to occur, we expect and then obtain:

$${[X_f,X_K]}_L=-X_{\{f,K\}}+\mathrm{div}\left(X_K\right)X_f.$$

(83)

Example 5. Robertson reactions

 

The reaction network:

$$X\stackrel{a}{\to }Y,Y+Y\stackrel{b}{\to }Y+Z,Y+Z\stackrel{c}{\to }X+Z$$

(84)

is a system of auto-catalytic reactions where $a,b,c$ are the reaction rates.

The mass action law gives a system clearly expressible in the gradient form ($\mathit{x}={(x,y,z)}^T$):

$$\dot{\mathit{x}}=B\left(\mathit{x}\right)\nabla H,B\left(\mathit{x}\right)=\left(\begin{array}{ccc}0\hfill & cyz+b{y}^2\hfill & -ax-b{y}^2\hfill \\ -cyz-b{y}^2\hfill & 0\hfill & ax\hfill \\ ax+b{y}^2\hfill & -ax\hfill & 0\hfill \end{array}\right)$$

(85)

with conserved $H=x+y+z$ (classical rule of mass conservation). Note additionally that $B\left(\mathit{x}\right)$ does not obey Jacobi’s identity, although its skew symmetry itself guarantees the conservation property (the system is almost Poisson) [12].

We are able to write down the system in a different form:

$$\dot{\mathit{x}}=\epsilon \nabla K,$$

(86)

ε being the totally anti-symmetric Cartesian-tensor of order three and:

$$K=-\frac{1}{2}a{x}^2-\frac{1}{3}b{y}^3-a{\int }_{y_0}^{y\left(t\right)}{x}^{\prime }d{y}^{\prime }-{\int }_{z_0}^{z\left(t\right)}(b{y}^{\prime 2}+c{y}^{\prime }{z}^{\prime })d{z}^{\prime };$$

(87)

therefore, we need a pair of reservoirs. In this form, the system is a Poisson one: the anti-symmetric structure matrix obeys the Jacobi identity; moreover, the system admits the Casimir function H, since it is obvious that $\epsilon \nabla H=0$; hence, we can proceed with the construction of the Darboux coordinates on a single symplectic leaf of the system, e.g., $(y,z)$, where $x=m_0-y-z$, $m_0$ constant.

The system reduces to:

$$\begin{array}{c}\dot{y}=\mu -ay-az-b{y}^2-cyz=-K_z,\hfill \\ \dot{z}=b{y}^2=K_y,\hfill \end{array}$$

(88)

where $\mu =a{m}_0$. To cast the above system in the gradient form, we need only a single reservoir:

$$K=-\mu z+\frac{1}{2}a{z}^2+\frac{1}{3}b{y}^3+{\int }_{z_0}^{z\left(t\right)}(a{y}^{\prime }+b{y}^{\prime 2}+c{y}^{\prime }{z}^{\prime })d{z}^{\prime }$$

(89)

and it is explicitly of the canonical form.

It is worth stressing that we can apply the Casimir function to the generator governing evolution of the system (87) to reduce the number of variables, but the formula will be different from that obtained applying the given Casimir to the equations of motion and then finding the reduced generator (89). The results are obviously unequal, but their differentials are cohomologically equivalent [1].

Example 6.

To illustrate that the ideas presented here work also in higher dimensions, we considered the five-dimensional system considered in [43].

$$\begin{array}{c}{\displaystyle \dot{p}=\alpha \left(z\right)-\beta py+\gamma pq,}\hfill \\ {\displaystyle \dot{q}=\delta \left(z\right)q-\epsilon qx-\gamma pq,}\hfill \\ {\displaystyle \dot{x}=\epsilon qx-\zeta x,}\hfill \\ {\displaystyle \dot{y}=\beta py-\eta y,}\hfill \\ {\displaystyle \dot{z}=-\alpha \left(z\right)p-\delta \left(z\right)q+\eta y+\zeta x.}\hfill \end{array}$$

(90)

Summing all the equations together, we find:

$$\frac{d}{dt}(p+q+x+y+z)=\alpha \left(z\right)(1-p),$$

(91)

and performing the Riemann–Stieltjes integral, we gain access to:

$$K=p+q+x+y+z-{\int }_{t_0}^t\alpha \left(z^{\prime }\right)(1-p^{\prime })d{t}^{\prime }=C_0$$

(92)

on the target level set.

Eliminating x, we obtain:

$$\begin{array}{c}{\displaystyle \dot{p}=\alpha \left(z\right)-\beta py+\gamma pq,}\hfill \\ {\displaystyle \dot{q}=\delta \left(z\right)q-\epsilon q\left(C_0+{\int }_{t_0}^t\alpha \left(z^{\prime }\right)(1-p^{\prime })d{t}^{\prime }-p-q-y-z\right)-\gamma pq,}\hfill \\ {\displaystyle \dot{y}=\beta py-\eta y,}\hfill \\ {\displaystyle \dot{z}=-\alpha \left(z\right)p-\delta \left(z\right)q+\eta y+\zeta \left(C_0+{\int }_{t_0}^t\alpha \left(z^{\prime }\right)(1-p^{\prime })d{t}^{\prime }-p-q-y-z\right).}\hfill \end{array}$$

(93)

which is writable in the Hamiltonian form with four reservoirs, provided we choose canonical pairs, for example $(p,y),(q,z)$.

Observe the interesting novelty: the traditional first integrals led to a reduction of the order of the ODE by one; effectively conserved quantities turned the problem into an integro-differential one instead. The significance of this fact and our capability to cope with it numerically will be discussed in future papers. The following example illustrates the possibility to shape the occurring systems according to our needs and taste.

Example 7.

Consider the pretty general system on ${\mathbb{R}}^3$:

$$\begin{array}{c}\dot{x}=f(x,y,z),\hfill \\ \dot{y}=g(x,y,z),\hfill \\ \dot{z}=h(x,y,z),\hfill \end{array}$$

(94)

with given initial conditions $(x_0,y_0,z_0)$.

If we would like to enclose the system as Poisson, we can assume:

$$\begin{array}{c}f=a{K}_y-b{K}_z,\hfill \\ g=-a{K}_x+c{K}_z,\hfill \\ h=b{K}_x-c{K}_y,\hfill \end{array}$$

(95)

with:

$$B\left(\mathit{x}\right)=\left(\begin{array}{ccc}0& a& -b\\ -a& 0& c\\ b& -c& 0\end{array}\right)$$

(96)

being the Poisson structure matrix, the entries of which satisfy the Jacobi identity:

$$a_z-b_y-a_x+c_z+b_x-c_y=0.$$

(97)

This constraint is the only requirement for $a,b,c$ to obey; despite this, they can be truly arbitrary.

Solving the linear system of Equation (95), we can determine the formal generator of motion up to a constant (of course, we demand that $K_x$, $K_y$, and $K_z$ be sufficiently smooth functions). Doing this, we have:

$$\begin{array}{c}{\displaystyle {\mathcal{K}}_x=\frac{1}{2}\left(\frac{h}{b}-\frac{g}{a}\right),}\hfill \\ {\displaystyle {\mathcal{K}}_y=\frac{1}{2}\left(\frac{f}{a}-\frac{h}{c}\right),}\hfill \\ {\displaystyle {\mathcal{K}}_z=\frac{1}{2}\left(\frac{g}{c}-\frac{f}{b}\right).}\hfill \end{array}$$

(98)

thus introducing the maximum three reservoirs (in fact, multiplicative factors could be chosen in a much more general way).

Obviously, for a given $B\left(\mathit{x}\right)$, we can introduce the formal Casimir function and then put the system in Hamiltonian form. Please notice that this general three-dimensional dynamical system is unique in a very specific sense: for dimensions greater than three, we obtain the whole family of non-potential generators equivalent in terms of gauge-type transformations. This remark leaves the door open for the future development of the presented reservoir approach.

6. Conclusions

Our article had a twofold objective. In Section 2 and Section 3, we reintroduced some notions of the geometric language of classical mechanics in convenient notation and summarized the basic terminology of kinetic chemistry and population dynamics. The considered examples were presented in a convenient Hamiltonian-like or Poisson-like framework.

In the main section of this paper, Section 4, we chose a rather pragmatic, not to say banausic, view of the theorem of the existence and uniqueness of the solutions of ordinary differential equations and exploited this for computational ends, to gain access to the so-called non-potential (or effective) integrals of motion. Here, we emphasize that they can be treated along the lines of a weak formulation of the integrals of motion (level manifolds). As a result, we proposed a new, para-Hamiltonian, description of ODEs.

Section 5 presented versatile examples showing how to use the constructed formalism. Among a variety of formal features, our approach is distinguished by the possibility of imposing the prescribed form of the problem under study, due only to our taste and convenience, which is clearly shown in Example 7.

Apparently, the innocent play with the dependent variables provided us with a tiny bit of additional information about the dynamical system not possessing first integrals in the classical sense. Non-potential, conserved quantities are of no use from the point of view of pure mathematics; yet, they can be very helpful in the numerical analysis of dynamical problems and the construction of new algorithms. In particular, we plan to construct geometric numerical integration schemes of any order for dissipative systems, analogous to those studied in [44]. The work in this direction is in progress.