# Schrödinger and Euler-Bernoulli beam equations: structure-preserving discretization and equivalence as port-Hamiltonian systems

Falak Zentout, Jesus-Pablo Toledo-Zucco, Lucie Baudouin, Denis Matignon

► To cite this version:

Falak Zentout, Jesus-Pablo Toledo-Zucco, Lucie Baudouin, Denis Matignon. Schrödinger and Euler-Bernoulli beam equations: structure-preserving discretization and equivalence as port-Hamiltonian systems. IEEE Conference on Decision and Control, Dec 2025, Rio de Janeiro (BRAZIL), Brazil. hal-05419201

HAL Id: hal-05419201

https://hal.science/hal-05419201v1

Submitted on 16 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Schrödinger and Euler-Bernoulli beam equations: structure-preserving discretization and equivalence as port-Hamiltonian systems

Falak Zentout

LAAS-CNRS,
Université de Toulouse, CNRS, INSA, UPS,
Toulouse, France.

falak.zentout@laas.fr

Lucie Baudouin

LAAS-CNRS,
Université de Toulouse, CNRS, INSA, UPS,
Toulouse, France.

lucie.baudouin@laas.fr

Jesus-Pablo Toledo-Zucco

LAAS-CNRS,
Université de Toulouse, CNRS, INSA, UPS,
Toulouse, France.

jesus-pablo.toledo-zucco@laas.fr

Denis Matignon

Fédération ENAC ISAE-SUPAERO ONERA,
Université de Toulouse,
Toulouse, France

denis.matignon@isae.fr

**Abstract**—In this work, the classical equivalence between Schrödinger and Euler-Bernoulli beam partial differential equations (PDEs) as closed physical systems is extended to open physical systems using the port-Hamiltonian framework: first, a damped version of both these models is presented; second, the possible collocated inputs and outputs of the two models are parameterized in the most general way; and third, a structure-preserving discretization method for a class of one-dimensional Boundary-Controlled Port-Hamiltonian System (BC-PHS) with collocated boundary actuation and sensing is developed. Unlike the classical Partitioned Finite Element Method (PFEM) approach, the proposed discretization method makes it possible to obtain an ordinary differential equation regardless of the boundary conditions, preventing the emergence of a differential-algebraic equation in the case of mixed boundary conditions. This novelty, recently validated for the wave and Timoshenko beam equations, is now extended to PDEs with a second-order differential operator. On the Schrödinger and Euler-Bernoulli equations, it is shown that the proposed numerical scheme can deal with a large type of boundary inputs and outputs.

**Index Terms**—Schrödinger, Euler-Bernoulli, port-Hamiltonian systems, partitioned finite element method

## I. INTRODUCTION

The Euler-Bernoulli beam model is a fundamental method for describing the bending behavior of beams and rods. It assumes that cross-sections remain flat and perpendicular to the beam's centerline during deformation while neglecting shear deformations, see e.g. [1]. Despite these simplifications, the model effectively captures the dynamic response of flexible

structures and is widely used in engineering applications, including bridges, aerospace structures, or other mechanical frameworks.

Similarly to the fundamental relation of dynamics in classical mechanics, the Schrödinger equation serves as a fundamental principle in quantum mechanics, governing the evolution of a system's wave function under a given potential. The solution of this equation is the density of probability of the massive non-relativistic particle under consideration to be in some given neighborhood. The applications of the Schrödinger equation are extensive, ranging from understanding the behavior of electrons in atoms [2] to the design of semiconductor devices. In chemical physics, it underpins the study of molecular structures and reaction rates [3]. Additionally, it plays a critical role in quantum computing research, where the manipulation of qubits relies on coherent evolution as governed by the Schrödinger equation [4].

The well-posedness of BC-PHSs in 1D has been shown through an algebraic condition on a rectangular matrix [5]. Recently, a necessary and sufficient algebraic condition for the well-posedness of second-order systems was established in [6] based on the inversibility of a boundary square matrix. Furthermore, asymptotic and/or exponential stability results can be achieved via output feedback even without internal dissipation [7]. In [8], it is established that the Schrödinger equation retains the generation of a contraction semigroup, along with the observability and admissibility properties, from the wave equation.

The equivalence between Schrödinger and Euler-
Bernoulli equations has been known for quite a long time, based on the property $(\partial_t - i\partial_{\xi\xi})(\partial_t + i\partial_{\xi\xi}^2) = \partial_{tt}^2 + \partial_{\xi\xi\xi\xi}^4$, see e.g. [9], and a renewed interest can be accounted for from the mathematical community, see [10], [11]. However the equivalence was only valid for closed physical systems, i.e. on $\mathbb{R}$, but could not take into account the possible interaction with an exterior domain through a control applied at the boundary.

In [12], the 1D Schrödinger equation is represented as a distributed Port-Hamiltonian System (PHS). The authors introduced a discretization that maintains the balance of the probability of presence equation at the discrete level; and since the computational domain is bounded, some appropriate equivalent boundary conditions have to be taken into account. In this work, simulations of Rabi oscillations to control the state of a system inside a quantum box are performed, and the numerical experiments include the transition between two levels of energy and the generation of Schrödinger cat states.

In the present paper, we adopt the port-Hamiltonian framework to emphasize boundary-controlled structure, energy preservation, and control-oriented formulation. The main contributions of this work are three-fold:

* We show that a unified formulation can deal with both the Schrödinger equation and the Euler-Bernoulli beam. This is well-suited to deal with existence and uniqueness of solutions and the design of static and dynamic control laws acting at the spatial boundaries of the domain.

* We introduce a damping model in both equations in an equivalent way and we are able to analyze it in the frequency domain.

* We develop a structure-preserving discretization method that preserves the energy behaviour, which is given in terms of the inputs and outputs as it is well-known for passive systems.

This paper is organized as follows. We detail the problem's setting in Section II, recalling the damped Schrödinger equation and the equivalent Euler-Bernoulli beam equation. In Section III, we recall the port-Hamiltonian representation of boundary controlled systems with second-order spatial derivative and we show the equivalences between the two equations. In Section IV, we propose a structure-preserving spatial discretization method and in Section V, we show some numerical simulations. Finally, in Section VI we give some conclusions and ideas of future works.

## II. PROBLEM STATEMENT

In this paper, we focus on analyzing the equivalence between two fundamental equations and their discretization using the PFEM within the PHS framework.

The damped Schrödinger equation under consideration is given by:

$$ \frac{\partial z}{\partial t}(\xi, t) = -i \frac{\partial^2 z}{\partial \xi^2}(\xi, t) - \alpha z(\xi, t) \quad (1) $$

where:

* $z(\xi, t)$ represents the complex-valued wave function (which is a density of probability),

* $\alpha \ge 0$ is a damping parameter, modeling constant potential loss effects.

On the other hand, the equivalent Euler-Bernoulli beam equation, modeling mechanical wave propagation in flexible structures, undergoes damping effects as well and a potential term, and takes the form:

$$ \frac{\partial^2 w}{\partial t^2}(\xi, t) = -\frac{\partial^4 w}{\partial \xi^4}(\xi, t) - 2\alpha \frac{\partial w}{\partial t}(\xi, t) - \alpha^2 w(\xi, t) \quad (2) $$

where:

* $w(\xi, t)$ is the transverse displacement of the beam,

* $\alpha \ge 0$ induces two combined effects: a damping of size $2\alpha$, and a shift of the eigenfrequencies of a coefficient $\alpha^2$.

The main goal of this work is to establish an equivalence between these two equations within the PHS formulation and to develop a structure-preserving numerical scheme for their discretization.

The discretization of infinite-dimensional systems often leads to numerical artifacts, such as energy drift or loss of stability. To overcome this, we use the PFEM, introduced in [13], which has shown to be a structure preserving approach. This approach ensures that the discrete model retains key properties of the continuous system, making it suitable for control and simulation.

## III. PORT-HAMILTONIAN FORMULATION AND EQUIVALENCES BETWEEN THE TWO EQUATIONS

In this section, we recall the port-Hamiltonian representation of a one-dimensional BC-PHS with boundary inputs and outputs. Then, we show that using above mentioned representation, one can either represent (1) and (2).

### A. BC-PHS with second order differential operator

We consider the following one-dimensional BC-PHS with internal damping, and with dynamical boundary inputs and outputs:

$$ \Sigma \begin{cases} \frac{\partial x}{\partial t}(\xi, t) = P_2 \frac{\partial^2 x}{\partial \xi^2}(\xi, t) - Gx(\xi, t), \\ \sigma(\xi, t) = \begin{bmatrix} x(\xi, t), \\ \partial_{\xi} x(\xi, t) \end{bmatrix}, \\ \begin{bmatrix} V_B \\ V_C \end{bmatrix} \begin{bmatrix} \sigma(b, t) \\ \sigma(a, t) \end{bmatrix} = \begin{bmatrix} U(t) \\ Y(t) \end{bmatrix} \end{cases} \quad (3) $$
with $\xi \in [a, b]$ being the spatial variable, $t > 0$ the time, $x(\xi, t) \in \mathbb{R}^2$ (for all $\xi \in [a, b]$ and $t > 0$) the state variable, and $\sigma(\xi, t) \in \mathbb{R}^4$ (for all $\xi \in [a, b]$ and $t > 0$) the co-energy variable and its first-order spatial derivative. The matrices $P_2$ and $G$ are defined as follows:

$$P_2 = \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix}, G = \begin{bmatrix} \alpha & 0 \\ 0 & \alpha \end{bmatrix}, \quad (4)$$

with $\alpha$ a positive scalar. Finally, the inputs and outputs $U(t) \in \mathbb{R}^4$ and $Y(t) \in \mathbb{R}^4$ are defined as a linear combination of the co-energy variables and its first-order spatial derivative evaluated at the spatial boundaries. The parametrization is achieved through the matrices $V_B \in \mathbb{R}^{4 \times 8}$ and $V_C \in \mathbb{R}^{4 \times 8}$, in such a way that (3) is passive with respect to the following storage energy or Hamiltonian:

$$H(t) = \frac{1}{2} \int_{a}^{b} x(\xi, t)^{\top} x(\xi, t) d\xi. \quad (5)$$

Indeed, one can show that the time derivative of the Hamiltonian is given by

$$\dot{H}(t) = \begin{bmatrix} \sigma(b,t) \\ \sigma(a,t) \end{bmatrix}^{\top} M \begin{bmatrix} \sigma(b,t) \\ \sigma(a,t) \end{bmatrix} - \int_{a}^{b} x^{\top}(\xi,t) G x(\xi,t) d\xi$$

with

$$M = \frac{1}{2} \begin{bmatrix} 0 & P_2 & 0 & 0 \\ P_2^{\top} & 0 & 0 & 0 \\ 0 & 0 & 0 & -P_2 \\ 0 & 0 & -P_2^{\top} & 0 \end{bmatrix}. \quad (6)$$

Thus, if we define the input $V_B$ and output $V_C$ matrices such that

$$\frac{1}{2} (V_B^\top V_C + (V_B^\top V_C)^\top) = M, \quad (7)$$

then one can verify the passivity of the Hamiltonian since

$$\dot{H}(t) = U(t)^{\top} Y(t) - \int_{a}^{b} x^{\top}(\xi, t) G x(\xi, t) d\xi. \quad (8)$$

**Remark 1:** System $\Sigma$ in (3) is a simplified version of a more general class of one-dimensional BC-PHSs in which spatial dependent parameters can be incorporated, namely mass densities and/or elasticity moduli.

## B. Equivalence with the Schrödinger Equation

System $\Sigma$ in (3) can be used to represent the Schrödinger Equation (1) with various types of boundary conditions. This can be achieved by decomposing the wave function $z(\xi, t)$ into its real and imaginary components and analyzing their dynamics. Indeed, if we define

$$x(\xi, t) = \begin{bmatrix} \Re\{z(\xi, t)\} \\ \Im\{z(\xi, t)\} \end{bmatrix},$$

one can show that the PDE in (3) is equivalent to (1). For simplicity, we use the notations $u(\xi, t) := \Re\{z(\xi, t)\}$

and $v(\xi, t) := \Im\{z(\xi, t)\}$. Since $z = u + iv$, then (1) becomes

$$\begin{aligned} \dot{u} + i\dot{v} &= -i(\partial_{\xi\xi}u + i\partial_{\xi\xi}v) - \alpha(u + iv), \\ &= -i\partial_{\xi\xi}u + \partial_{\xi\xi}v - \alpha u - i\alpha v, \\ &= \partial_{\xi\xi}v - \alpha u + i(-\partial_{\xi\xi}u - \alpha v), \end{aligned}$$

which gives us the following representation

$$\frac{\partial}{\partial t} \begin{bmatrix} u \\ v \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} \frac{\partial^2}{\partial \xi^2} \begin{bmatrix} u \\ v \end{bmatrix} - \begin{bmatrix} \alpha & 0 \\ 0 & \alpha \end{bmatrix} \begin{bmatrix} u \\ v \end{bmatrix} \quad (9)$$

corresponding to the PDE of system $\Sigma$ in (3).

1) Energy Analysis: The Hamiltonian of the Schrödinger equation (1) is given by:

$$H_z(t) = \frac{1}{2} \int_a^b z(\xi, t) \bar{z}(\xi, t) d\xi.$$

Its time derivative, using integrations by parts and the Schrödinger equation in (1), gives:

$$\begin{aligned} \dot{H}_z(t) &= \Re \left( \int_a^b z_t(\xi, t) \bar{z}(\xi, t) d\xi \right) \\ &= \Re \left( \int_a^b (-i \partial_{\xi\xi} z(\xi, t) - \alpha z(\xi, t)) \bar{z}(\xi, t) d\xi \right) \\ &= \Re \left[ -i \partial_{\xi\xi} z(\xi, t) \bar{z}(\xi, t) \right]_a^b - 2\alpha H_z(t) \\ &= -\frac{i}{2} \left[ \partial_{\xi\xi} z(\xi, t) \bar{z}(\xi, t) - \overline{\partial_{\xi\xi} z(\xi, t) z(\xi, t)} \right]_a^b - 2\alpha H_z(t) \end{aligned}$$

so that $\dot{H}_z(t) = \begin{bmatrix} z(b,t) \\ \partial_{\xi\xi}z(b,t) \\ z(a,t) \\ \overline{\partial_{\xi\xi}z(a,t)} \end{bmatrix}^* \mathbb{M} \begin{bmatrix} z(b,t) \\ \partial_{\xi\xi}z(b,t) \\ z(a,t) \\ \overline{\partial_{\xi\xi}z(a,t)} \end{bmatrix} - 2\alpha H_z(t)$,

with: $\mathbb{M} := \frac{1}{2} \begin{bmatrix} 0 & P_2 & 0 & 0 \\ P_2^* & 0 & 0 & 0 \\ 0 & 0 & 0 & -P_2 \\ 0 & 0 & -P_2^* & 0 \end{bmatrix}$ and $P_2 = -\frac{i}{2}$.

In this latter case, the input and output matrices $V_B$ and $V_C$ such that $V_B^* V_C = \mathbb{M}$ with $V_B, V_C \in \mathbb{C}^{2 \times 4}$.

## C. Equivalence with the Euler-Bernoulli beam

Similarly to Section III-B, we show now that (3) can be used to represent the Euler-Bernoulli beam equation (2) with different types of boundary conditions. This can be achieved by defining the following state

$$x(\xi, t) = \begin{bmatrix} \partial_{\xi\xi} w(\xi, t) \\ \partial_t w(\xi, t) + \alpha w(\xi, t) \end{bmatrix}.$$

Taking the time derivative of the previous state, we get

$$\begin{aligned} \dot{x} &= \begin{bmatrix} \partial_{\xi\xi} \dot{w} \\ \dot{w} + \alpha \dot{w} \end{bmatrix} = \begin{bmatrix} \partial_{\xi\xi} (\dot{w} + \alpha w) - \alpha \partial_{\xi\xi} w \\ -\partial_{\xi\xi} \xi \xi w - 2\alpha \dot{w} - \alpha^2 w + \alpha \dot{w} \end{bmatrix} \\ &= \begin{bmatrix} \partial_{\xi\xi} (\dot{w} + \alpha w) - \alpha (\partial_{\xi\xi} w) \\ -\partial_{\xi\xi} \xi \xi w - \alpha (\dot{w} + \alpha w) \end{bmatrix}. \end{aligned}$$

Thus, if we define here $u := \partial_{\xi\xi}w$ and $v := \partial_t w + \alpha w$, then we obtain (9) again, which is equivalent to the PDE of $\Sigma$ in (3).
IV. STRUCTURE-PRESERVING SPATIAL
DISCRETIZATION

In this section, we propose a structure-preserving
spatial discretization method that preserves the energy
structure of the initial PDE in (3). This means that
the Hamiltonian (5) and its balance equation in (8)
are preserved at the discrete level. Firstly, we write
down the weak form of (3). Secondly, we assign
a port-Hamiltonian structure to the finite-dimensional
approximation. Finally, we show that the Hamiltonian
of the discretized model mimics the Hamiltonian of
the initial PDE. The work follows [13, §. 6.3], where
second-order differential operators are involved, but
takes advantage of a general parametrization of the
boundary conditions, as done originally in [14] for
first-order differential operators, with the advantage
of ensuring an explicit ODE in the end, instead of
an implicit Differential-Algebraic Equation (DAE): this
approach is strongly related to the Hellinger-Reissner
method, which was already analyzed in [15].

A. Weak form on the discretized space

We consider the Galerkin finite element method. The
state variables $x_1(\xi, t)$ and $x_2(\xi, t)$ of the PDE in
(3) are approximated using the interpolation functions
$\phi_1 \in H^2([a, b], \mathbb{R}^{2N_1})$ and $\phi_2 \in H^2([a, b], \mathbb{R}^{2N_2})$,
respectively in the following way:

$$
\begin{bmatrix} x_1(\xi, t) \\ x_2(\xi, t) \end{bmatrix} \approx \begin{bmatrix} x_1^{ap}(\xi, t) \\ x_2^{ap}(\xi, t) \end{bmatrix} := \begin{bmatrix} \phi_1(\xi)^{\top} x_1^d(t) \\ \phi_2(\xi)^{\top} x_2^d(t) \end{bmatrix}, \quad (10)
$$

where $x_i^{ap}$ is the approximation of $x_i$, and $x_i^d \in \mathbb{R}^{2N_i}$
(for $i = \{1, 2\}$) is the time-dependent discretized state
vector that contains the discrete values of $x_i$ and its
corresponding spatial derivative $\partial_\xi x_i$. We replace the
approximation (10) in the PDE in (3)

$$
\begin{align*}
\phi_1(\xi)^{\top} \dot{x}_1^d(t) &= \partial_{\xi\xi} \phi_2(\xi)^{\top} x_2^d(t) - \alpha \phi_1(\xi)^{\top} x_1^d(t), \\
\phi_2(\xi)^{\top} \dot{x}_2^d(t) &= -\partial_{\xi\xi} \phi_1(\xi)^{\top} x_1^d(t) - \alpha \phi_2(\xi)^{\top} x_2^d(t).
\end{align*}
$$

In the following, we avoid spatial and time dependence
since it is clear from the context. Multiplying, on the
left, the first equation by $\phi_1$ and the second by $\phi_2$, and
integrating in $\xi \in [a, b]$, one gets

$$
\begin{bmatrix} E_1 & 0 \\ 0 & E_2 \end{bmatrix} \begin{bmatrix} \dot{x}_1^d \\ \dot{x}_2^d \end{bmatrix} = \begin{bmatrix} -R_1 & D_{12} \\ D_{21} & -R_2 \end{bmatrix} \begin{bmatrix} x_1^d \\ x_2^d \end{bmatrix} \quad (11)
$$

with

$$
E_i = \int_a^b \phi_i \phi_i^\top d\xi, \quad R_i = \int_a^b \phi_i \alpha \phi_i^\top d\xi, \quad i = \{1, 2\}, \\
D_{12} = \int_a^b \phi_1 \partial_{\xi\xi} \phi_2^\top d\xi, \quad D_{21} = -\int_a^b \phi_2 \partial_{\xi\xi} \phi_1^\top d\xi.
$$

B. Port-Hamiltonian representation of the finite-
dimensional system

We define the matrices $\Phi(\xi) = \begin{bmatrix} \phi_1(\xi) & 0 \\ 0 & \phi_2(\xi) \end{bmatrix}$ and $\Psi = [\Phi^{(b)} \ \partial_{\xi} \Phi^{(b)} \ \Phi^{(a)} \ \partial_{\xi} \Phi^{(a)}]$, of respective sizes $(2N_1 + 2N_2) \times 2$ and $(2N_1+2N_2) \times 8$. Then, applying integration by parts to $D_{12}$ and $D_{21}$ twice, one can show that

$$
\begin{bmatrix} 0 & D_{12} \\ D_{21} & 0 \end{bmatrix} = 2\Psi M \Psi^{\top} + \begin{bmatrix} 0 & -D_{12}^{\top} \\ -D_{12}^{\top} & 0 \end{bmatrix} \quad (12)
$$

with $M$ defined in (6). Now, we approximate the inputs
and outputs defined in (3) as follows:

$$
U = V_B \begin{bmatrix} x^{(b)} \\ \frac{\partial_{\xi} x^{(b)}}{x^{(a)}} \\ \frac{\partial_{\xi} x^{(a)}}{x^{(a)}} \end{bmatrix} \approx V_C \begin{bmatrix} \Phi^{(b)} x^{(a)} \\ \frac{\partial_{\xi} \Phi^{(b)} x^{(a)}}{\Phi^{(a)} x^{(a)}} \\ \frac{\partial_{\xi} \Phi^{(a)} x^{(a)}}{\Phi^{(a)} x^{(a)}} \end{bmatrix} = V_B \Psi^{\top} x^{(a)},
$$

$$
Y = V_C \begin{bmatrix} x^{(b)} \\ \frac{\partial_{\xi} x^{(b)}}{x^{(a)}} \\ \frac{\partial_{\xi} x^{(a)}}{x^{(a)}} \end{bmatrix} \approx V_C \begin{bmatrix} \Phi^{(b)} x^{(a)} \\ \frac{\partial_{\xi} \Phi^{(b)} x^{(a)}}{\Phi^{(a)} x^{(a)}} \\ \frac{\partial_{\xi} \Phi^{(a)} x^{(a)}}{\Phi^{(a)} x^{(a)}} \end{bmatrix} = V_C \Psi^{\top} x^{(a)}.
$$

If we define the input matrix as

$B := \Psi V_C^\top,$

then, denoting

$$
E = \begin{bmatrix} E_1 & 0 \\ 0 & E_2 \end{bmatrix}, \quad M_d = \begin{bmatrix} 0 & D_{12} \\ D_{21} & 0 \end{bmatrix}, \quad R = \begin{bmatrix} R_1 & 0 \\ 0 & R_2 \end{bmatrix},
$$

the discretized model (11) can be written as

$$
\begin{align*}
E \dot{x}^d &= (M_d - R)x^d + BU - BU \\
&= (M_d - R)x^d + BU - B V_B \Psi^\top x^d \\
&= (M_d - R)x^d + BU - \Psi V_C^\top V_B \Psi^\top x^d \\
&= (M_d - \Psi V_C^\top V_B \Psi^\top - R)x^d + BU.
\end{align*}
$$

Finally, we define the interconnection matrix as

$$
J := M_d - \Psi V_C^\top V_B \Psi^\top \quad (13)
$$

and we show that it is a skew-symmetric matrix ($J = -J^\top$). To this end, we prove $J + J^\top = 0$. Notice that from (12), one obtains $M_d + M_d^\top = 2\Psi M \Psi^\top$. Then, using (13), (7), and (12), one can show:

$$
\begin{align*}
J + J^{\top} &= M_d - \Psi V_C^{\top} V_B \Psi^{\top} + M_d^{\top} - \Psi V_B^{\top} V_C \Psi^{\top}, \\
&= M_d + M_d^{\top} - \Psi (V_C^{\top} V_B + V_B^{\top} V_C) \Psi^{\top}, \\
&= M_d + M_d^{\top} - 2\Psi M \Psi^{\top} = 0.
\end{align*}
$$

Finally, the discretized version of the BC-PHS in (3) is
given by the following PHS :

$$
\Sigma_d \begin{cases} E \dot{x}^d(t) = (J - R)x^d(t) + BU(t), \\ Y(t) = B^\top x^d(t). \end{cases} \quad (14)
$$

in which $J = -J^\top$ and $R = R^\top \ge 0$.

Remark 2: Notice that one can also consider the
internal damping parameter as a spatial dependent
function as soon as $\alpha(\xi) \ge 0$ for all $\xi \in [a, b]$. When $\alpha$
is constant, the matrix $R$ becomes $R = \alpha E$.
### C. Energy structure of the discretized model

We define the energy of the discretized model in (14) as the Hamiltonian in (5) replacing $x$ by $x^{ap}$ from (10)

$$ 
\begin{aligned}
H_d(t) &:= \frac{1}{2} \int_a^b x^{ap}(\xi, t)^\top x^{ap}(\xi, t) d\xi, \\
&= \frac{1}{2} \int_a^b x^d(t)^\top \Phi(\xi) \Phi(\xi)^\top x^d(t) d\xi = \frac{1}{2} x^d(t)^\top E x^d(t).
\end{aligned}
$$

Now, computing its energy balance through the dynamic (14), we get

$$ 
\begin{aligned}
\dot{H}_d &= \frac{1}{2} (x^{d\top} E \dot{x}^d + \dot{x}^{d\top} E x^d) \\
&= x^{d\top} E \dot{x}^d = x^{d\top} ((J - R)x^d + BU) \\
&= -x^{d\top} Rx^d + x^{d\top} BU = -x^{d\top} Rx^d + Y^{\top} U,
\end{aligned}
$$

in which we recover the same structure of the initial energy balance in (8).

### D. Hermite Polynomials

Following [16, §. 6.2], since the variables at the boundaries can involve the first spatial derivatives, we use the Hermite polynomials $\varphi_i$ and $\psi_i$ for the test function. In this case, we use the following polynomials of degree 3:

$$ 
\begin{aligned}
\varphi_i(\xi) &= \alpha_i + \beta_i \xi + \gamma_i \xi^2 + \delta_i \xi^3, \\
\psi_i(\xi) &= a_i + b_i \xi + c_i \xi^2 + d_i \xi^3,
\end{aligned}
$$

with $i = \{a, b\}$ with the following characteristics:

$$ 
\begin{aligned}
\varphi_a(a) &= 1, & \varphi_b(a) &= 0, & \psi_a(a) &= 0, & \psi_b(a) &= 0, \\
\varphi_a(b) &= 0, & \varphi_b(b) &= 1, & \psi_a(b) &= 0, & \psi_b(b) &= 0, \\
\varphi'_a(a) &= 0, & \varphi'_b(a) &= 0, & \psi'_a(a) &= 1, & \psi'_b(a) &= 0, \\
\varphi'_a(b) &= 0, & \varphi'_b(b) &= 0, & \psi'_a(b) &= 0, & \psi'_b(b) &= 1,
\end{aligned}
$$

as illustrated in Fig. 1.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Left: Hermite polynomials of a single element. Right: First derivative of the Hermite polynomials.</figcaption></figure>

## V. NUMERICAL SIMULATION

In this section, we show the numerical simulation of a cantilever Euler-Bernoulli beam and the probability density of a Schrödinger equation. In the first case, the corresponding boundary conditions are mixed, whereas in the second one, the boundary conditions are of the type Dirichlet. In both cases we consider $N = 100$

test functions for each state and for each first spatial derivative of the states. This leads to an ODE of 400 states. The Matlab codes and the corresponding animations are provided in a GitLab project<fnref n="1" />.

### A. Euler-Bernoulli beam: Cantilever beam

We consider the Euler-Bernoulli beam in (2) with unitary length, $[a, b] = [0, 1]$, with low internal dissipation $\alpha = 0.1$ and the following boundary conditions:

$$ 
\begin{aligned}
\partial_t w(a, t) &= 0, & \partial_\xi \partial_t w(a, t) &= 0, \\
\partial_\xi \xi w(b, t) &= -k_1 \partial_\xi \partial_t w(b, t) + \tau_{ext}(t), \\
\partial_\xi \xi \xi w(b, t) &= k_2 \partial_t w(b, t) + F_{ext}(t),
\end{aligned}
$$

where the first two equations correspond to the left clamped side of the beam and the other two equations correspond to an external force $F_{ext}(t)$ and external torque $\tau_{ext}(t)$ which are applied to the right side. The parameters $k_1 = 0.2$ and $k_2 = 0.2$ are the boundary damping coefficients related to the angular and transversal velocity, respectively.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Simulation of a cantilever beam</figcaption></figure>

**Simulation:** The objective is to move the beam from the zero equilibrium position to a new equilibrium position. The initial conditions are set to zero, i.e., $w(\xi, 0) = 0$ and $\partial_t w(\xi, 0) = 0$. The external force and torque are the following functions $F_{ext}(t) = 0.4(1 - e^{-5t})$ and $\tau_{ext}(t) = 0.2(1 - e^{-5t})$. In Fig. 2, we show the temporal response of the beam. The initial condition in violet arrives to the final equilibrium in orange. We remark that the exact equilibrium of this example when $\alpha = 0$ is $w^*(\xi) = (\tau^* - F^*b) \frac{\xi^2}{2} + F^* \frac{\xi^3}{6}$ with $\tau^* = 0.2$, and $F^* = 0.4$, which is close to the one shown in orange in Fig. 2. The oscillations during the transient are due to the values of $k_1$ and $k_2$. By modifying these values, the oscillations can be reduced. The animation can be found in the Gitlab project<fnref n="2" />

<footnote n="1">¹https://gitlab.com/ToledoZucco/discretization_EB_SE.</footnote>

<footnote n="2">²https://gitlab.com/ToledoZucco/discretization_EB_SE/-/blob/main/Beam_Animation.gif</footnote>
## B. One-dimensional Schrödinger equation

The simulation of the Schrödinger equation for $\alpha = 0$ shows the evolution of $|z(\xi, t)|^2$, which gives the probability of finding a particle at the point $\xi$, at time $t$. To this end, we use homogeneous Dirichlet boundary conditions (fixed ends), i.e., $z(a, t) = 0$ and $z(b, t) = 0$.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Simulation of the Schrödinger equation.</figcaption></figure>

**Simulation:** The objective is to study the evolution of the probability density for a given initial condition $\Re\{z(\xi, 0)\} = \frac{1}{(\sigma\sqrt{2\pi})}e^{-(\xi-0.5)^2/(2\sigma^2)}$ and $\Im\{z(\xi, 0)\} = 0$ with $\sigma = \sqrt{0.005}$. In Fig. 3, we can see the evolution of the probability density. One can see that at the initial condition (in violet), the most probable point in the spatial domain is the $\xi = 0.5(m)$. Then at $t \approx 0.04(s)$, there is a very low probability at $\xi = 0.5(m)$ and a high probability at $\xi = 0.25(m)$ or $\xi = 0.75(m)$. Finally, at $t \approx 0.08(s)$, the probability density approaches the initial condition (in orange) and since $\dot{H}(t) = 0$, the curve repeats the same profile in time. The animation can be found in the Gitlab project<fnref n="3" />.

## VI. CONCLUSIONS AND FUTURE WORK

A significant aspect of this work is the employment of the PFEM for structure-preserving discretization, which guarantees that the numerical approximation upholds the Hamiltonian properties of the original continuous equations. We validated this discretization method through numerical simulations. Moreover, the equivalence between the two PDEs has been extended to boundary-controlled systems, relating the collocated boundary controls and observations.

While these findings are promising, they pave the way for further exploration. Several challenges remain, providing exciting avenues for future research: one promising research direction is extending the port-Hamiltonian approach to non-collocated inputs and outputs. Addressing this requires alternative formulations that preserve the system's energy-based structure. One potential solution is the Scattering Passive Systems framework, recently explored in infinite-dimensional systems [17].

## ACKNOWLEDGEMENTS

This work has been partially funded by the French National Research Agency (ANR), via project IMPACTS<fnref n="4" />, grant number ANR-21-CE48-0018.

## REFERENCES

[1] A. Brugnoli, D. Alazard, V. Pommier-Budinger, and D. Matignon, “Port-Hamiltonian formulation and symplectic discretization of plate models Part II: Kirchhoff model for thin plates,” *Applied Mathematical Modelling*, vol. 75, pp. 961–981, 2019.

[2] D. J. Griffiths and D. F. Schroeter, *Introduction to Quantum Mechanics*. Cambridge University Press, Aug. 2018.

[3] I. N. Levine, *Quantum chemistry*. Boston: Pearson, 2014.

[4] M. A. Nielsen and I. L. Chuang, *Quantum Computation and Quantum Information: 10th Anniversary Edition*. Cambridge University Press, Dec. 2010.

[5] Y. Le Gorrec, H. Zwart, and B. Maschke, “Dirac Structures and Boundary Control Systems Associated with Skew-Symmetric Differential Operators,” *SIAM Journal on Control and Optimization*, vol. 44, pp. 1864–1892, 2005.

[6] B. Elghazi, B. Jacob, and H. Zwart, “Well-posedness of a class of infinite-dimensional port-Hamiltonian systems with boundary control and observation,” in *Proc. 5th IFAC Workshop on Control of Systems Governed by Partial Differential Equations (CPDE)*, IFAC, IEEE-CSS, June, 2025. https://arxiv.org/abs/2501.16862.

[7] B. Augner and B. Jacob, “Stability and stabilization of infinite-dimensional linear port-Hamiltonian systems,” *Evolution Equations and Control Theory*, vol. 3, no. 2, pp. 207–229, 2014.

[8] H. Zwart, Y. Le Gorrec, and B. Maschke, “Relating systems properties of the wave and the Schrödinger equation,” *Evolution equations and control theory*, vol. 4, no. 2, pp. 233–240, 2015.

[9] J.-D. Polack, “Viscous interpretation and time formulation of hysteretic damping,” *Acta Acustica united with Acustica*, vol. 87, no. 3, pp. 352–358, 2001.

[10] I. V. Volovich, “On the equivalence between the Schrödinger equation in quantum mechanics and the Euler-Bernoulli equation in elasticity theory,” *p-Adic Numbers, Ultrametric Analysis and Applications*, vol. 17, no. 1, pp. 78–84, 2025.

[11] S. N. Gavrilov, A. M. Krivtsov, and E. V. Shishkina, “Energy transport in a free Euler-Bernoulli beam in terms of Schrödinger’s wave function,” *Mechanics Research Communications*, p. 104382, 2025.

[12] G. Verrier, G. Haine, and D. Matignon, “Modelling and structure-preserving discretization of the Schrödinger as a port-Hamiltonian system, and simulation of a controlled quantum box,” in *International Conference on Geometric Science of Information*, pp. 392–401, Springer, 2023.

[13] F. L. Cardoso-Ribeiro, D. Matignon, and L. Lefèvre, “A partitioned finite element method for power-preserving discretization of open systems of conservation laws,” *IMA Journal of Mathematical Control and Information*, vol. 38, no. 2, pp. 493–533, 2021.

[14] J.-P. Toledo-Zucco, D. Matignon, C. Poussot-Vassal, and Y. Le Gorrec, “Structure-preserving discretization and model order reduction of boundary-controlled 1D port-Hamiltonian systems,” *Systems & Control Letters*, vol. 194, p. 105947, 2024.

[15] A. Brugnoli, G. Haine, and D. Matignon, “Explicit structure-preserving discretization of port-Hamiltonian systems with mixed boundary control,” *IFAC-PapersOnLine*, vol. 55, no. 30, pp. 418–423, 2022.

[16] G. Allaire and A. Craig, *Numerical Analysis and Optimization: An Introduction to Mathematical Modelling and Numerical Simulation*. Oxford University Press, 2007.

[17] J.-P. Toledo-Zucco, D. Matignon, and C. Poussot-Vassal, “Scattering-passive structure-preserving finite element method for the boundary controlled transport equation with a moving mesh,” *IFAC-PapersOnLine*, vol. 58, no. 6, pp. 292–297, 2024.

<footnote n="3">³https://gitlab.com/ToledoZucco/discretization_EB_SE/-/blob/main/Schrodinger_Animation.gif</footnote>

<footnote n="4">⁴https://impacts.ens2m.fr/</footnote>