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

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 ω 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: given on some open domain with ; obviously, an overset dot stands for a shortcoming denoting time derivation.

Given , the existence and continuity of the solutions are guaranteed. The assumption of 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 such that: where denotes the mentioned integral. Note that B is not determined uniquely, since we can add any solution of the homogeneous equation to the particular solution:The exterior product could be defined as: being explicitly anti-symmetric (note that in , the structure provided by such a product is equivalent to the one given by the usual cross-product). The Formula (3) explicitly demands . However, if we are already given the problem in the form of (2), the assumption of 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 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 , we have:(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: where: and the Lie algebra of functions defined on the phase space.

This also gives rise to the first-order differential operator: known as a Hamiltonian vector field.

Now, we call it the the first integral generator of motion if: understood componentwise. In other words: or simply: denoting by ⌋ the substitution of a vector field into the form (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: 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 , we have .

It surely obeys the continuity law: or, in other words:

Now, because the phase fluid is incompressible: on a target energy level set. Now, let us ponder , and then, we vary the constant as to obtain the condition: where c is a rightful constant.

The fluid should also undergo Euler’s equation (the fluid’s particle velocity denoted by ; here, we considered just the fluid of an ensemble of free harmonic oscillators with the equations of motion ): which in an elementary manner leads to Bernoulli’s law, if we take the inner product of both sides with and perform trivial integration: from which we can evaluate the pressure function and constant appearing in the formula, as the pressure needs to be non-negative.

3. -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:

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;

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: would undergo conditions 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: with a properly adjusted set of constants (note that the above equations fulfill the previously mentioned restrictions, if 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: where denotes the number of prey caught by the predator per unit of time and 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 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 -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 as the basic object (hence the name -systems). It serves as a generator for the equations of motion:

For instance, in the case of the Lotka-Volterra system (20) we have: with the generator: 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]).

Similarly, we can formulate the Hamiltonian-like

The set of chemical reactions is termed the reaction network. Having: 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 ( 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 , 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]): 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 , 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 in all their arguments in a given domain, in a classical sense.

Given a Pfaff form where the path of the integration is along this trajectory, denoted by 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). where and can be derived by the comparison of both equations. Hence:In other words, we define in the phase space one-form , which is neither closed nor exact:

There is ongoing energy exchange with the environment (not necessarily a loss) due to non-potential function . 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 , and theoretically, it is a redundant variable; hence: 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”: also giving rise to one-form , 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.

The quantity: is conserved.

We may easily check: 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 ). Since, for general , we have , 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:

In the future, we will denote by the expression standing next to in .

Therefore, we can define: 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, becomes the Pfaffian form).

Since (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: . Motion begins with ; energy as a number at a given time is invariant with respect to the change of coordinates (“frozen” time transformations); since the initial value is also invariant, it induces the invariance of the integral while changing the description of the system so that: where the summation over repeating indices is assumed.

Demanding the para-Hamiltonian framework to hold in both sets of symplectic coordinates: we obtain the condition for action differentials to differ by an exact differential (as usual): where the time-dependent differentials cancel due to the common initial conditions. Traditionally performing proper Legendre transforms, we can express 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 ; the left-over condition gives: where we tacitly assumed . 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 to x or from to 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: giving rise to the Poisson bracket: 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: moreover: which is indeed the case.

We can regard the condition: for any 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:

Adjoining external forces to the language in which we describe the system, we obtain a closed system; hence, , and hence, the non-potential force exerts a kind of “pressure” on the ensemble of systems in the phase space (see below);

If , then we view and as equal, so reaching for the phase pressure obtained below for the phase fluid and considering phase volume as , we can pick the monotonically growing component of the expression 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 such that , but , 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: seeing that compressible terms ( is of course a tangent vector to the phase flow, namely ) 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 , 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: as a symplectic form. During the former considerations, it was preserved by the flow of Hamiltonian vector field , and now, it satisfies: as it provided by the equations of motion: , hence , so multiplying both sides by p, we obtain .

is defined unambiguously throughout the phase space as a section over , 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 . Hence, , 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: 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: yielding: and, accordingly:

Since we are playing with ODE system, when we assume that , 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: where the reservoirs are defined by: so that .

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: hence: 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 , so C may depend on the K value, so on the phase trajectory of interest. Here, as always, .

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

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

Since all fluid particles obey these equations of motion, the continuity equation yields: where c is a constant.

Writing Bernoulli’s law: and remembering that , we see:

Fortunately, we know solutions to the damped oscillator, being ; hence, . Thus, we that see there is no danger of the variable part of pressure growing to infinity. Provided that the engaged constants obey: where 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.

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 [

The non-potential generator of motion is: where the reservoir variable is given by: turning K into an effectively conserved quantity.

Note that taking the 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

Note that onward, to avoid cumbersome writing in integrals such as (

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

We claim that substrates

A quick look at these confirms there is no resource function where the reservoir variables are given by: so that

Note that again, in the which is preserved due to the courtesy of K being constant during the evolution.

Note that (69) has a fixed point at . This equilibrium is unstable when ; if , it is stable. The case 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 , then the Hartman–Grobman theorem on linearization often fails at predicting the existence of the center [19,20,21].

We put the thermostatted system in the form: which gives: 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

The almost-Poisson structure provided in [ 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 (

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.

The Lotka–Volterra model describes the basic predator–prey interaction (with a linear response): 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: claiming that

We should observe that such a formulated LV problem is given on the phase space and, as a Poisson system, has a compressible phase fluid: ; 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 , where the equations of motion become: with the incompressible phase fluid on the symplectic phase space with the separable resource function:

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

Moreover, we have: 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.

The transition from the Poisson dynamics governed by the resource function

This remark is easily verifiable on a case-by-case basis, provided that the Poisson bracket of with coordinates is preserved.

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

The reaction network: is a system of auto-catalytic reactions where

The mass action law gives a system clearly expressible in the gradient form ( with conserved

We are able to write down the system in a different form: ε being the totally anti-symmetric Cartesian-tensor of order three and: 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

The system reduces to: where 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].

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

Summing all the equations together, we find: and performing the Riemann–Stieltjes integral, we gain access to: on the target level set.

Eliminating x, we obtain: which is writable in the Hamiltonian form with four reservoirs, provided we choose canonical pairs, for example

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.

Consider the pretty general system on with given initial conditions

If we would like to enclose the system as Poisson, we can assume: with: being the Poisson structure matrix, the entries of which satisfy the Jacobi identity: This constraint is the only requirement for

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

Obviously, for a given

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.