Energetic-Property-Preserving Numerical Schemes for Coupled Natural Systems

: In this paper, we propose a method for deriving energetic-property-preserving numerical schemes for coupled systems of two given natural systems. We consider the case where the two systems are interconnected by the action–reaction law. Although the derived schemes are based on the discrete gradient method, in the case under consideration, the equation of motion is not of the usual form represented by using the skew-symmetric matrix. Hence, the energetic-property-preserving schemes cannot be obtained by straightforwardly using the discrete gradient method. We show numerical results for two coupled systems as examples; the ﬁrst system is a combination of the wave equation and the elastic equation, and the second is of the mass–spring system and the elastic equation.


Introduction
Recently, against the background of complication of designing industrial products, the necessity of coupled simulation of plural phenomena according to different equations has increased. To carry out such simulations, first, it is necessary to construct and discretize appropriate mathematical models before numerical calculations. However, even if the equations for individual phenomena are known, naive combination of the equations easily results in violating important law of physics such as the energy-conservation law [1,2].
In this paper, we propose energetic-property-preserving numerical schemes for coupled models of two given natural systems. Popular models for coupled systems are the Port-Hamiltonian systems [3]. From the viewpoint of geometric mechanics, in this approach, the equations are described by using the Dirac structure, which is a geometric structure generalized from the symplectic structure and the Poisson structure [4]. This structure describes the energy transfer between the systems, thereby deriving the model equation in such a way that the energetic property of the coupled systems is preserved. It is also known that the resultant coupled system is also port-Hamiltonian and hence models for large coupled systems can be derived in a systematic way. In addition, the Lagrangian formalism of this approach is also reported [5][6][7].
Regarding numerical schemes for the port-Hamiltonian systems, energy-preserving and passivity-consistent numerical schemes were proposed by Celledoni and Høiseth [8]. Focusing on a discrete energy balance equation and the stability under interconnection, Celledoni and Høiseth proposed a definition of discrete port-Hamiltonian systems. Port-Hamiltonian systems are typically written as dx dt = (B(x) − R(x))∇H(x) + G(x)u, x(0) = x 0 , y = G(x) ∇H(x), (1) where x ∈ R n is the state variables, u ∈ R m is the input, and y ∈ R m is the output. B(x) is a skew-symmetric matrix, which is often corresponding to the Hamiltonian structure of the underlying systems. G(x) is a n × m matrix. H : R n → R is typically an energy function. R(x) is a positive definite matrix, which is corresponding to the dissipation property of the system. For conservative systems, R(x) = 0. In [8], two conservative systems, one of which is supposed to be given while the other is a controller to be designed, are assumed to be interconnected by u =ỹ,ũ = −y (2) so that the total energy is preserved: In fact, the above system is shown to be d dt with a skew-symmetric matrix C. By using this structure, energy-preserving and passivity-consistent numerical schemes are successfully derived.
In this paper, we consider a slightly different situation where two natural systems are interconnected by the action-reaction law. Natural systems are systems that are derived from Lagrangian defined as the difference between the kinetic energy and the potential energy. This type of coupled systems naturally arise in industrial simulations. The systems under consideration can be infinite dimensional and hence we consider three types of interconnections, that is, those between an ODE and an ODE, a PDE and a PDE, and an ODE and a PDE. As explained in Section 2, due to the action-reaction law, the systems are not always written in the form of Equation (4). Hence, a different approach to designing energetic-property-preserving numerical method is required for these systems, which we propose in this paper. This paper is organized as follows. In Section 2, we explain the coupled systems and the energy conservation law of these systems. Then, we propose energetic-property-preserving numerical schemes for the coupled systems. Section 3 deals with the systems consisting of finite dimensional systems and Section 4 does the systems including infinite ones. Numerical examples are shown in Section 5.

Remark 1.
The above three types of coupled systems often appear in industrial simulations. For example, let us consider simulations of the stringed musical instruments, such as the guitar and the piano [9]. The one side of strings of these instruments are typically attached to a bar-shaped part of the instruments, which is called the bridge. Each string and the bar are modeled as infinite-dimensional natural systems and they are in contact at a single point. Hence, this case corresponds to the interaction between two PDEs. Meanwhile, the strings of the piano are excited by the hammer, which is modeled as the interaction between a wave-type PDE (a model of each string) and a mass-spring ODE (a model for the hammer) (see Figure 1). In addition, designing numerical schemes for the above two cases often reduces to the case of interconnections between two ODEs by semi-discretization in space (cf. Section 4).

Figure 1.
An example of coupled systems (a rough model of a piano). In this example, the string and the bridge are modeled as PDEs, which they are in contact at a single point. Excitation of the string by a hammer is typically modeled as an ODE.

Interconnection of Natural Systems
In this paper, we consider interconnection of two natural systems with dissipation terms. Natural systems are systems derived from the standard Lagrangian that is defined as the difference between the kinetic energy and the potential energy. For example, where m is a constant, u is the vector of the state variables, and V is the potential function, is a natural system because this is the Euler-Lagrange equation associated with the Lagrangian The system in Equation (5) often arises with the damping term in which G is a positive semidefinite operator. If G is proportional to the identity operator, the operator G essentially denotes the damping coefficient. For simplicity of notation, regarding G as a 2-form, we often write for vectors or functions w 1 and w 2 G(w 1 , w 2 ) = w 1 Gw 2 (8) or G(w 1 , w 2 ) = w 1 Gw 2 dx.
These are nonnegative if w 1 = w 2 .
Another example of a natural system is the Euler-Lagrange partial differential equation of the following form: where ρ is typically the density and δH/δu the variational derivative. The variational derivative is the gradient in functional spaces, which is defined by the function that satisfies dH is the Fréchet derivative and ·, · L 2 is the standard L 2 inner product. In fact, the partial differential Equation (10) has the Lagrangian Similar to Equation (5), Equation (10) is often accompanied by dissipation terms: Example 1. For example, the wave equation on the interval (0, L) under the periodic boundary condition is the Euler-Lagrange partial differential equation, of which Lagrangian is In this example, H is and the variational derivative of this functional is because under the periodic boundary condition it follows that from the formal calculation Throughout this paper, we suppose that two natural systems with dissipation terms of the form in Equation (7) or Equation (13) are in contact at a single point and through this point the two systems are exerting a force on each other. As noted in the Introduction, we consider the following three cases: interconnections between an ODE and an ODE, a PDE and a PDE, and an ODE and a PDE.
To begin with, we consider the first case. As mentioned above, we assume that systems are interconnected in such a way that the force between each two subsystems is described by the action-reaction law. More specifically, we assume that the coupled model is written in the form u 1 ∈ R n 1 , u 2 ∈ R n 2 are vectors of state variables which are possibly of different dimensions. e j is the jth unit vector representing the contact point. f is the scalar function of t, which we define so that the total energy of this system is preserved when the dissipation terms do not exist. First, the time derivative of the total energy of the system is d dt Substitution of the equation yields where ( du k dt ) i denotes the ith component of the vector du k dt . In the absence of the dissipation terms, this becomes Thus, for conservation of the energy is sufficient. For this condition to be held, we impose that also for the systems with dissipation terms. From Equation (20) and the first equality of Equation (25), we get which gives As a result, the coupled model equation is Remark 2. If we write Equation (28) without the dissipation terms can be written in the Hamiltonian formalism: where B is the matrix that has only 4 non-zero components Because, unfortunately, the matrix B is not skew-symmetric, the equation is not of the form in Equation (4); however, as shown in the next theorem, the total energy is conserved.
for any solutions to the equations of the interconnected systems the energy is not increasing: In particular, in the absence of the dissipation terms, that is, G 1 = G 2 = 0, the energy is preserved: Proof. From ith and jth elements of Equation (33), the following equations hold: Substitution of Equation (34) into the above yields which is equivalent to In the same way, we have and, hence, Thus, we get for all t. Hence, using the assumption on the initial conditions in Equation (32), for any t, we have Meanwhile, as shown in Equation (22), the time derivative of the energy is computed as below: Therefore, from Equation (38), we obtain In the same way, we can obtain similar theorems for the other two types of the interconnected systems. Firstly, we consider where f is a scalar function of t. δ is the Dirac delta function, and hence the equalities in the above equations are in the weak sense; however, the following computations are naively performed for simplicity. The terms with the Delta function represent that the two systems are interconnected at x = x 1 for the first system and x = x 2 for the second. For this coupled system, the total energy is defined by Theorem 2. Suppose that the boundary condition is appropriately imposed so that Provided that for any solutions to the equations of the interconnected systems the energy is not increasing: In particular, in the absence of the dissipation terms, the energy is preserved: Secondly, we consider of which the total energy is Theorem 3. Suppose that the boundary condition is appropriately given so that Provided that for any solutions to the equations of the interconnected systems the energy is not increasing: In particular, in the absence of the dissipation terms, the energy is preserved: We omit the proofs of Theorems 2 and 3 because they are proved exactly in the same way as Theorem 1.

Energetic-Property-Preserving Numerical Schemes for the Finite Dimensional Coupled Systems
In this section, we propose energetic-property-preserving numerical schemes for the above finite-dimensional coupled systems. Although the equation is not of the form in Equation (4), we use the discrete gradient similarly to CelledoniHøiseth [8]. A discrete gradient is defined as follows (see, e.g., [10]) The first condition is a discrete counterpart of the chain rule The second condition assures that a discrete gradient is indeed an approximation of the gradient. A discrete gradient is not uniquely determined and hence several ways for designing a discrete gradient have been proposed. A typical way is the average vector field method [11]: Another choice for obtaining a discrete gradient is the automatic discrete differentiation, which is a method to automatically derive a discrete gradient in a similar way to the automatic differentiation; see [12].
In what follows, we use the following symbols for the discrete systems. The step sizes are denoted, for example, by ∆t and ∆x. An approximation of u(n∆t) is denoted by u (n) , where u (n) can be a vector Similarly, for partial differential equations, an approximation of u(n∆t, j∆x) is denoted by u (n) j . For simplicity of notation, we also use the following difference operators: and the averaging operator Note that the difference operators and the averaging operator commute; for example, We first consider the coupled system in Equation (33). For the other two types of systems, the numerical schemes are obtained by replacing the discrete gradient by the discrete variational derivative [1,13] or by discretizing the partial differential equation in space in such a way that the resultant semi-discrete scheme becomes a finite dimensional natural system; we discuss this below.
We discretize the Hamiltonian form of Equation (33) in the following way: Under an assumption that corresponds to Equation (32), this numerical scheme preserves the energetic properties of the systems in the following sense.
for any solutions to the numerical scheme in Equation (61), the total energy is not increasing: In particular, in the absence of the dissipation terms, that is, G 1 = G 2 = 0, the energy is preserved: Proof. The proof is almost the same as the continuous case. First, from the definition of f (n+ 1 2 ) , it follows that and hence In the same way, we get Therefore, we obtain for all n from the assumption on the initial conditions in Equation (62). Meanwhile, the time difference of the discrete energy function is Substitution of the numerical scheme in Equation (61) yields

Energetic-Property-Preserving Numerical Schemes for the Coupled Systems Incorporating the Infinite Dimensional Systems
The numerical scheme for the partial differential Equation (40) can be obtained in the following way. First, the equation is semi-discretized in space into the ordinary differential equation of the form in Equation (20), which can be achieved in the straightforward way by discretizing the potential function in space as seen in the following example. For a general method of semi-discretization, see, e.g., [1,11,13].

Remark 3.
For example, the wave equation with the wave speed c and the damping coefficient γ on the interval (0, L) under the periodic boundary condition is discretized to where u j (t) is an approximation to u(t, j∆x), ∆x = L/N. This is an ordinary differential equation of the form Second, we need to discretize the term with the Dirac delta function: Because the delta function is defined by requiring for any continuous function g, we need to discretize this functional while satisfying this condition in some sense. A straightforward discretization would be where j 0 is the integer that has the nearest j 0 ∆x to x 0 . A more sophisticated way is first approximating the Dirac function by a continuous function, such as the Gauss function and then discretizing this function; for example, with a small σ. c 0 is the normalization constant that is defined by requiring ∑ j ε j ∆x = 1.
Following the procedure described above, we can semi-discretize the system of the partial differential equations in Equation (40) into and Equation (47) into where ε, ε 1 , ε 2 are vectors that approximate the Dirac functions andˆover the symbols denotes approximated values by the semi-discretization. The following theorem is obtained in the same way as shown in Section 2.
Theorem 5. Provided that for any solutions to the semi-discretized scheme the energy is not increasing: where µ 1 , µ 2 are 1 or ∆x, which are determined depending on whether the system under consideration originates from an ODE or a semi-discretization of a PDE. In the absence of the dissipation terms, the energy is conserved. Proof. In the same way as before, we have Hence, the proof is competed if we have which follows from the assumption in Equation (76) on the initial conditions and We show this equality in what follows. Computing the innerproduct of µ 1 ε 1 and the equations forû 1 gives Substitution off yields and hence we get In the same way, we also have forû 2 and which shows Equation (80).
Theorem 6. Provided that for any solutions to the numerical scheme the total energy is not increasing: In particular, in the absence of the dissipation terms, that is,Ĝ 1 =Ĝ 2 = 0, the energy is preserved: Proof. This theorem is obtained in the same way as before. The time difference of the discrete energy function is Substitution of the numerical scheme in Equation (82) yields Thus, the proof is completed if we can show that which follows from under the assumption in Equation (81) on the initial conditions.
From the scheme in Equation (82) forû (n) 1 , we have which shows In the same way, we have and hence This shows Equation (86).

Remark 5.
The accuracy of the proposed schemes depends on the choice of the discrete gradient and also on the spatial discretization if the systems are infinite dimensional. If the discrete gradient is symmetric, then the scheme is second-order accurate in temporal direction because the scheme is also symmetric; otherwise, the scheme is first-order accurate. For example, the average vector field method in Equation (56) yields a symmetric discrete gradient. As for the spatial direction, the accuracy mainly depends on the discretization of the Dirac delta function; however, since the equations are discretized by using the finite difference method, it is actually difficult to determine the accuracy of the discretization of this functional because the Dirac delta function does not admit the standard Taylor expansion. For precise analysis, the finite element method should be applied.

Remark 6.
Because the above schemes are implicit, they are stable but computationally expensive. There are several techniques to reduce the computational effort for solving the system of equations, which arises from the implicit scheme. For example, the predictor-corrector method proposed in [1] deserves consideration.

Numerical Example
In this section, we show two numerical examples. The first example is a coupled system that consists of the wave equation and the elastic equation under the boundary conditions The total energy of this system is Following Equation (61), we discretize these equations into where we use for simplicity the same step size ∆x and ∆t for both systems and k 1 and k 2 correspond to the interconnecting points: k 1 ∆x = x 1 , k 2 ∆x = x 2 . (u 1 ) (n) and (u 2 ) (n) are the vectors that consist of (u 1 ) (n) j s and (u 2 ) (n) j s. D 2 and D 4 are the difference matrix defined by which are the approximations of ∂ 2 u ∂x 2 , ∂ 4 u ∂x 4 by the central difference operators under the boundary conditions in Equation (90). Without the dissipation terms, each of these discrete systems is the Hamiltonian equation of which energy function is, respectively, where N 1 and N 2 are the numbers of the nodes. We set L 1 = L 2 = 1, ∆x = 1/100, ∆t = 0.001. The initial conditions are set as The positions of the interconnection are x 1 = x 2 = 0.2. We first set the damping coefficients γ 1 and γ 2 to 0. The numerical solution from t = 0 to t = 1.0 are shown in Figures 2 and 3. We observe that the wave packet shown in the initial state of the wave equation splits and propagates towards the connecting point and the right boundary. Then, at the connecting point, the wave packet is reflected while a slow wave is excited in the elastic equation.
Although this causes an energy exchange, the total energy is preserved up to the rounding errors, as shown in Figure 4. In addition, in this example, if the energy is bounded, then (v 1 ) (n) L 2 are all bounded. In fact, the stability of the numerical scheme follows from these conservation laws because combined with the Dirichlet boundary condition and the boundedness of δ + x (u 1 ) (n) L 2 and δ + x (u 2 ) (n) L 2 shows that (u 1 ) (n) L ∞ and (u 2 ) (n) L ∞ are also bounded from the discrete Poincaré inequality [13]. Second, we set both the damping coefficients γ 1 and γ 2 to 10. The numerical solution and the total energy from t = 0 to t = 1.0 are shown in Figures 5, 6, and 7, respectively. In particular, from the results in Figure 7, it is confirmed that the total energy is monotonically dissipated.   (88) with the damping terms. The excitation due to the energy exchange is smaller than that in the case without the damping terms. The second example is a combination of the mass-spring system and the elastic equation:

Concluding Remarks
In this paper, we design energetic-property-preserving numerical methods for coupled natural systems. For such methods for coupled systems, Celledoni and Høiseth proposed energy-preserving numerical methods [8] aiming at controlling port-Hamiltonian systems. In this paper, we consider similar but slightly different models; restricting the target systems to natural systems, we coupled two systems by using the law of action-reaction. When the two systems with no dissipation term are coupled, it is essentially required that the time derivative of the state variables of the two systems are equal at the connecting point for conservation of the total energy. The model and the energetic-property-preserving numerical methods are constructed by determining the interaction between the systems so that this condition holds. As shown in the first numerical experiment, the conservation of energy often leads to the stability of numerical schemes. Such a systematic approach for deriving stable schemes is useful for simulations of large coupled systems because it is generally difficult to derive stable numerical methods for the complex systems.
For future work, the proposed technique should be applied for deriving stable numerical schemes for complicated systems that actually appear in industrial applications. In addition, because the proposed approach in this paper is based on the discrete gradient method, this approach derives implicit methods, which are not often feasible for large-scale systems. Therefore, it is important to develop a method for relaxing the derived implicit schemes to be explicit. For symplectic integrators, Tao proposed a method for rewriting any Hamiltonian system to a separable Hamilton system so that explicit numerical schemes are derived by using the splitting method [14]. It is desirable to develop explicit energetic-property-preserving numerical methods for coupled systems by using such a method. Acknowledgments: The authors are grateful to the anonymous referees for their useful suggestions and comments that were helpful for improving the paper.

Conflicts of Interest:
The authors declare no conflict of interest.