Fereshteh Asgari
email: fereshteh.asgari@telecom-sudparis.eu

Alexis Sultan
email: alexis.sultan@telecom-sudparis.eu

Haoyi Xiong

Mounim Vincent Gauthier
email: vincent.gauthier@telecom-sudparis

El Yacoubi

Mounim El-Yacoubi
email: mounim.el_yacoubi@telecom-sudparis

CT-Mapper: mapping sparse multimodal cellular trajectories using a multilayer transportation network

Keywords: Mobile phone, Mobile networks signaling, Multimodal transportation network, HMM, Unsupervised learning, Mobile trajectories mapping, Intelligent transportation systems

come

algorithm that enables the mapping of mobile phone traces over a multimodal transport network. One of the main strengths of CT-Mapper is its capability to map noisy sparse cellular multimodal trajectories over a multilayer transportation network where the layers have different physical properties and not only to map trajectories associated with a single layer. Such a network is modeled by a large multilayer graph in which the nodes correspond to metro/train stations or road intersections and edges correspond to connections between them. The mapping problem is modeled by an unsupervised HMM where the observations correspond to sparse user mobile trajectories and the hidden states to the multilayer graph nodes. The HMM is unsupervised as the transition and emission probabilities are inferred using respectively the physical transportation properties and the spatial coverage of antenna base stations. To evaluate CT-Mapper we collected cellular traces with their corresponding GPS trajectories for a group

Introduction

Macroscopic analysis of the traffic flow in large metropolitan areas is a challenging task. This is especially true when multiple transit authorities are in charge of different transport networks (road, train, subway). Due to the lack of a common source of information across these transit systems, it is often hard for city authorities to grasp a unified view of mobility patterns. In this context, mobile phone data have recently become an attractive source of information about mobility behavior. Thanks to the ubiquitous usage of mobile phones, mining mobile phone data becomes a promising way to understand multimodal human mobility [START_REF] Reddy | Using mobile phones to determine transportation modes[END_REF][START_REF] Doyle | Utilising mobile phone billing records for travel mode discovery[END_REF][START_REF] Smoreda | Spatiotemporal data from mobile phones for personal mobility assessment[END_REF] ranging from identifying a mobile user daily path to recording transportation usage (e.g., taking train, metro, bus, etc.) in a large metropolitan area. Traditional approaches of mobility studies used GPS to accurately sense spatial data with a localization error bound ≤ 50m. Although they ensure the collection of fine-grained mobility trajectories (as shown in Fig. 1b), GPS-based data collection has two main drawbacks: first, it causes high energy consumption, and second, it is constrained to a limited group of users (e.g. taxi drivers [START_REF] Kang | Exploring human movements in singapore: A comparative analysis based on mobile phone and taxicab usages[END_REF] or a group of car drivers [START_REF] Giannotti | Unveiling the complexity of human mobility by querying and mining massive trajectory data[END_REF]). GPS sensing, therefore, is not suitable for collecting large-scale data from metropolitan area populations. By contrast, cellular data provided by network operators does not suffer from these issues, and has become recently, as a result, a new source of mobility information. Signaling information from mobile network operators (CDRs -Call Data Records-) has been used as a valuable source of mobility information for large scale population [START_REF] Smoreda | Spatiotemporal data from mobile phones for personal mobility assessment[END_REF][START_REF]Detecting mobility patterns in mobile phone data from the ivory coast[END_REF][START_REF] Agarwal | Large scale model for information dissemination with device to device communication using call details records[END_REF].

Localization of mobile phone users with antennas (i.e., cellular towers), nonetheless, provides only coarse-grained mobility trajectories at antenna level, with a varying localization error of hundred meters in densely populated cities, and within several kilometers in rural areas [START_REF] Smoreda | Spatiotemporal data from mobile phones for personal mobility assessment[END_REF]. Given the resulting cellular mobility trajectories (i.e., a sequence of antenna id s) and the location of each antenna as shown in Fig. 1c, it might be difficult to observe the road or metro station that the user passes by (as shown in Fig. 1a).

In order to collect cellular mobility trajectories using mobile phones, previous works [START_REF]Detecting mobility patterns in mobile phone data from the ivory coast[END_REF][START_REF] Smoreda | Spatiotemporal data from mobile phones for personal mobility assessment[END_REF] usually extracted the trajectories from Call Detailed Records (CDR), where the CDR of a user restores the antenna id and the time-stamp of each of his/her mobile calls. To understand human mobility, these works were mostly limited to aggregating the trajectories from a user's long-term CDR data in order to determine the frequently-visiting locations and the visiting time (e.g., the park he/she usually passes by during the 07:00-09:00 window of working days). As such, the techniques proposed by previous works are not suitable for estimating the precise mobility trajectories on the road/transportation network using the CDR cellular trajectories. Furthermore, one sample of CDR data (i.e., one call record) can be obtained only when the user places a call, making human mobility data between two consecutive calls irretrievable, especially when the time duration between the two calls is long (e.g., the inter-call mobility between the two calls in Fig. 1d).

Thus, even though it has been studied widely, CDR is unlikely to be a good data source for the trajectory mapping problem. Considering the time sparsity drawbacks of CDRs, we use, in this work, a new passive capturing technique to efficiently extract the position of the base stations the mobile phone is connected to. This technique analyzes the signaling channel of the data mobile network in order to extract base station locations. This way of capturing the mobility of users is salable and provides a higher sampling rate than CDR-based sensing.

The sparse cellular trajectories are collected and provided upon the request of the experiment participants to the network operator. Considering privacy issues [START_REF] De Montjoye | Unique in the crowd: The privacy bounds of human mobility[END_REF], the network operator localizes each mobile user using an antenna id, and further records each user's antenna id with time-stamp periodically (e.g., every 15 minutes in our study). Compared to the user's real trip (in Fig. 1a), the sparse cellular trajectory (in Fig. 1e) partially measures the user's mobility with coarse-grained localization. The objective of our work is to map each sparse cellular trajectory1 into the multimodal transportation network, in order to obtain the sequence of network nodes that the user passes by. For example, given the cellular trajectory shown in Fig. 1e and the transportation network of the Paris metropolitan area shown in Fig. 2, our goal is to recover the sequence of nodes of the real trip shown in Fig 1a.

The common approach for mapping cellular trajectories into the metropolitan transportation (usually road) network is to first collect a large amount of cellular trajectories and then to manually label each cellular trajectory with the corresponding intersection sequence, an intersection being a graph node associated with a junction between two roads. The next phase is to train a supervised mobility model (e.g., HMM) using the labeled cellular trajectories, in order to build a probabilistic model mapping antenna id sequences to intersection sequences. After training, given a new user cellular trajectory, the supervised model predicts, as the mapping result, a sequence of intersections, having the maximal likelihood of generating the antenna id sequence. However, obtaining the labeled cellular trajectories to cover the road/transportation networks and all cellular towers of a metropolis is not practical, as it costs too much human efforts for trajectory collection and labeling. We propose, in this paper, to solve the cellular trajectory mapping problem using an unsupervised mobility model, that does not require collecting and labeling any trajectories.

Given the above examples and target research goals, the key issues in designing the unsupervised mobility model include:

1) Given the antenna id sequence in a cellular trajectory, retrieve the se-quence of road/rail intersections that the user passes by given a database storing the multimodal transportation network -The transportation network covering and connecting multiple types of transportation modes (e.g., rail, metro, highway , etc.) is named multimodal transportation network [START_REF] Liu | Data Model and Algorithms for Multimodal Route Planning with Transportation Networks[END_REF], in which each node is either a road intersection or a station of a rail transportation mode (i.e., subway, tramway and train), and each edge is a connection between intersections (e.g., the pathway connecting a metro station and a bus stop). Obviously, it is nontrivial to extract the precise user path from the multimodal transportation network using the antenna id sequences. In summary, the main contributions of this work are:

• We propose to study the problem of mapping cellular trajectories to the multimodal transportation network, in order to obtain the precise mobility of the users. To the best of our knowledge, this is the first work addressing these issues. In particular, rather than mapping cellular trajectories using the supervised mapping algorithms with labeled mobility data, we propose to use an unsupervised mapping algorithm leveraging on topological properties of the transportation network, so as to eliminate the tedious human labeling efforts in building the mobility model. • We collect real cellular trajectories of a group of users in the Paris metropolitan area with the help of a French telecom operator, then evaluate our mapping algorithm using the data. Through the extensive evaluation with cellular trajectories covering more than 2500 intersection nodes and 3 physical layers , 1000 metro and subway stations, we show that our algorithm maps the cellular trajectory onto the multimodal transportation network of the Paris metropolitan area and with good accuracy given the sparsity of user cellular trajectories. This algorithm also achieves up to 20% higher accuracy compared to a baseline approach, that exploits for unsupervised HMM parameter estimation, the complexity and topology of the multilayer network, without considering the transportation properties of network edges.

The rest of this paper is structured as follows: Sec.2 presents related work.

Sec. 3 gives an overview of the proposed system. Sec. 4 presents the details of the unsupervised estimation of HMM parameters and explains how the two main probability distributions used for mapping are derived. In Sec. 5 , we evaluate our proposed algorithm and the paper ends by a discussion and a conclusion in Sec. 6.

Related Work

General Human Mobility Models

A considerable amount of Human Mobility studies have been devoted to analyze trajectories of individuals based on their traces. Spatial characteristics such as center of the mass, radius of gyration and statistical characteristics revealed a number of scaling properties in human trajectories: Gonzalez et al in [START_REF] González | Understanding individual human mobility patterns[END_REF] and Brockmann et al [START_REF] Brockmann | The scaling laws of human travel[END_REF] showed a truncated power-law tendency in the distribution of jump length. It was observed that most individuals travel only over a short distance, and there is only a few who travel regularly over hundred kilometers. Further studies [START_REF] Simini | A universal model for mobility and migration patterns[END_REF][START_REF] González | Understanding individual human mobility patterns[END_REF] showed that travel patterns collapse into a single spatial probability distribution, indicating that, despite the diversity of their travel history, humans follow simple reproducible patterns. In addition, statistical analysis confirms that individuals' movement follows spatio-temporal patterns [START_REF] Giannotti | Unveiling the complexity of human mobility by querying and mining massive trajectory data[END_REF][START_REF] Giannotti | Trajectory pattern mining[END_REF][START_REF] Csáji | Exploring the mobility of mobile phone users[END_REF] which can help defining mobility models. In all mentioned studies, multimodal mobility aspects were not taken into account. One objective, in this work, is to investigate the mobility patterns of trajectories through the multimodal transportation network and to explore how these patterns are affected by the multiplicity of the layers of the network. Early mobility studies relied on expensive data collection methods, such as surveys and direct observation. Trajectories were mostly defined as Origin-Destination (OD), and were mapped over the desirable graph to retrieve an optimal path solution which is usually the shortest path between the Origin and Destination [START_REF] Simini | A universal model for mobility and migration patterns[END_REF][START_REF] Yuan | Discovering regions of different functions in a city using human mobility and pois[END_REF][START_REF] Giannotti | Unveiling the complexity of human mobility by querying and mining massive trajectory data[END_REF][START_REF] Giannotti | Trajectory pattern mining[END_REF]. Although recent studies have been trying to infer the traffic flow using additional traffic data [START_REF] Abadi | Traffic flow prediction for road transportation networks with limited traffic data[END_REF], they still fail to retrieve the real path taken by individuals.

Mapping Algorithms

Along with mobility studies, applications such as navigation systems, traffic monitoring and public transportation tracking, used GPS data to track individuals or any moving object [START_REF] Thiagarajan | Accurate, low-energy trajectory mapping for mobile devices[END_REF][START_REF] Goh | Online map-matching based on hidden markov model for real-time traffic sensing applications[END_REF][START_REF] Newson | Hidden markov map matching through noise and sparseness[END_REF][START_REF] Hu | Adaptive kalman filtering for vehicle navigation[END_REF][START_REF] Xu | Development and application of a kalman filter and gps error correction approach for improved map matching[END_REF][START_REF] Hunter | Scaling the mobile millennium system in the cloud[END_REF]. A variety of statistical approaches such as Expectation Maximization (EM) [START_REF] Hunter | Scaling the mobile millennium system in the cloud[END_REF], Kalman Filter [START_REF] Hu | Adaptive kalman filtering for vehicle navigation[END_REF][START_REF] Xu | Development and application of a kalman filter and gps error correction approach for improved map matching[END_REF] and Hidden Markov Model (HMM) [START_REF] Thiagarajan | Vtrack: Accurate, energy-aware road traffic delay estimation using mobile phones[END_REF][START_REF] Thiagarajan | Accurate, low-energy trajectory mapping for mobile devices[END_REF][START_REF] Goh | Online map-matching based on hidden markov model for real-time traffic sensing applications[END_REF][START_REF] Newson | Hidden markov map matching through noise and sparseness[END_REF][START_REF] Hummel | Map matching for vehicle guidances[END_REF] were used to map noisy sequential location data over transportation networks. Most of these mapping algorithms have used GPS data as they provide accurate location data with an error of about 50 meters. Moreover, using labeled data, supervised models were considered and trained to optimize model parameters in an automatic way. Once the models are trained, they are used to find the most likely path in the network assigned to sequences of noisy location data. However, most of these mapping algorithms were developed to map noisy data over road networks without considering other mobility modes.

Human Mobility Modeling with CDR Cellular Trajectories

Recently, because of the expeditious growth of mobile phones, Call Data Records (CDR) have been providing great data sets for human mobility studies as they are collected continuously for all active cellular phones. CDRs, however, have two significant limitations: first, they are sparse in time because they are generated only when a phone engages in a voice call or text message exchange;

and second, they are coarse in space and less precise than GPS location data, because they record location only at the granularity of a cellular antenna (with an average error of 175 meter in dense populated areas and up to 2 kilometers in less denser areas). Nonetheless, the fact that almost the entire population is already equipped with cell phones [START_REF] Smoreda | Spatiotemporal data from mobile phones for personal mobility assessment[END_REF] allows for studying important aspects of individual mobility such as inferring transportation modes. Cellular network data were, for instance, used to classify different transportation modes for longdistance travels [START_REF] Smoreda | Spatiotemporal data from mobile phones for personal mobility assessment[END_REF][START_REF] Doyle | Utilising mobile phone billing records for travel mode discovery[END_REF]. Thiagaran et al. in [START_REF] Thiagarajan | Accurate, low-energy trajectory mapping for mobile devices[END_REF] used cellular signal data with combination of cellphone sensors to develop a supervised mapping algorithm in order to overcome the limitation of GPS data. While previous works have used cellular data to map long trajectories, this work proposes an unsupervised mapping algorithm that maps the sparse cellular trajectories over the multimodal transportation network in the Paris metropolitan area. This approach could be used for large scale smart-phone users for further studies in traffic estimation.

Such a mapping is important for the development of smart cities and smart mobility.

Studies of smart cities in the past were limited to analyzing multimodal transportation networks without considering large scale real mobility data. The main goal of multimodal mobility studies is to improve public transportation monitoring and to reduce traffic congestion [START_REF] Liu | Data Model and Algorithms for Multimodal Route Planning with Transportation Networks[END_REF][START_REF] Van Nes | Design of multimodal transport networks : a hierarchical approach[END_REF][START_REF] Aguiar | Leveraging electronic ticketing to provide personalized navigation in a public transport network, Intelligent Transportation Systems[END_REF] . Considering the aforementioned observations and the fact that the majority of trajectory mapping problems are developed for mono-modal transportation networks (specifically road networks), we believe that there is a gap in the literature. This study aims at bridging this gap by mapping cellular sparse data of smartphones over the multimodal transportation network in the Ile-de-France metropolitan area.

The multimodal mapping results may help not only optimizing the multimodal transportation network, but also investigating the multimodal mobility behavior of individuals in metropolitan areas.

CT-Mapper System Overview

In this section, we first formulate the search problem of CT-Mapper, and introduce the datasets collected for mapping. We then analyze the computational complexity of the mapping problem over the collected datasets, and finally present the framework of CT-Mapper.

Problem Statement

In this section, we first formulate the problem by defining several key concepts used in our approach. Function Ψ indicates the layer of each node Ψ :

Definition 1. Multilayer Transportation Graph

V → L in G. Transportation Layer G l = (V l , E l ) is a subset of G, where V l = {v|v ∈ V, Ψ(v) = l} and E l = {< v i , v j >∈ E, Ψ(v i ) = Ψ(v j ) = l}. Each node v i
is characterized by its latitude and longitude (i.e., the geographical position CrossLayer edge set E cl ⊂ E defines the edges with pair of nodes not belonging to the same layer:

v i =< lat, lon > i )
E cl = {< v i , v j >∈ G|Ψ(v i ) ̸ = Ψ(v j )}
The multilayer Transportation graph is characterized by its adjacency matrix issue in GTP protocol can be found in [START_REF] Xu | Acculoc: Practical localization of performance measurements in 3g networks[END_REF]. Finally, when the mobile looses the signal or is turned off, a message closing the session is sent. With modern smartphone applications that emit and receive data on a regular basis (i.e.

W ij ∈ R |V |×|V | .
email, push notification), it is expected that the GTP tunnel for a given user remains constantly maintained, enabling us to sample the user position at each network event (handover and radio resource allocation).

GPS Trajectory Data -

To evaluate the accuracy of our proposed mapping algorithm, GPS data were used as ground truth. A group of participants were asked to install the "Moves" smartphone application [START_REF]Moves[END_REF] to record their GPS locations. The GPS locations provided by "Moves" were analyzed to extract real trajectories of participants. The underlying transportation network used in this study is the multimodal transportation network of Ile-de-France which is modeled by several separated graph layers corresponding each to a different transportation mode, interconnected together into a multiplex G. To build this multiplex, multiple geospatial datasets, namely the road network from the National Geographic Institute (IGN)[29] and the rail transport network (train and metro) from Open-StreetMap (OSM) [START_REF]Openstreetmap project[END_REF] were aggregated. Each node in G is either a road intersection, a rail station or a metro station. A key feature of the proposed multimodal transportation network is its modeling of transitions between different transport modes during a given trip. Cross-layer transition modeling are ensured by adding CrossLayer appropriate edges between layers.

Computational Complexity of the

Although such a multilayer representation of the transportation network enables us to model and define trajectories using different transportation modes, it also increases the complexity of the underlying network. To highlight this fact we use the "search complexity " metric to show how difficult is it to find the sequence of segments that compose the truth path over the map. This metric describe how hard is it to find this sequence of segments by chance. if the sequence is difficult find by chance then the same apply for an algorithm also (regardless of its internal design).

First in Table 1 we illustrates some topological differences between each layer in the multilayer graph G. For example, the average length between two consecutive intersections is rather heterogeneous across different transportation layers. To quantitatively assess network complexity, we use entropy measure to characterize the ease/difficulty of navigation in a network using "the search information" developed in [START_REF] Sneppen | Hide-and-seek on complex networks[END_REF], and in [START_REF] Rosvall | Networks and cities: An information perspective[END_REF]. In order to summarize the work of [START_REF] Sneppen | Hide-and-seek on complex networks[END_REF] and [START_REF] Rosvall | Networks and cities: An information perspective[END_REF] the search information entropy is the Shannon entropy of the probability of finding a given destination (in Eq. 1) by chance, higher is the entropy more difficult it will be for any search algorithm to find a right destination.

Eq. ( 1) defines the probability for a random walker starting at node s with degree k s to reach node t. Consequently, in Eq. ( 2) we define the search entropy of a graph as the sum over all shortest paths {SP st } from node s to node t in G averaged over all possible pair of nodes (s, t) in graph G. As a result, by computing the average entropy of all the possible paths in G, we can express the relative complexity (S avg ) of finding a given path in a given graph G.

P [SP st ] = 1 k s ∏ j∈SPst 1 k j -1 (1)
S avg = 1 N (N -1) N ∑ s=1 N ∑ t=1 -log 2 ∑ {SPst} P [SP st ] (2)
In Fig. 4a we plot the average entropy of each layer of the multimodal transportation graph of Ile-de-France along with the average entropy of the interconnected multilayer network. We observe that the average entropy is As it shows clearly, the complexity of the multilayer graph is higher than each of its layers taken separately, regardless of its size. We define σ = S/ log 2 (N ) as the average graph path's entropy relative to its size and δ = (S avg -S R )/ log 2 (N ) to describe how a graph compares with its random counterpart in terms of its node degree, irrespective of the network size.

As a conclusion, the search complexity of finding the right path in the multilayer transportation graph increases compared with a single layer graph. This is due to two effects: firstly when different layers are combined together in a multilayer graph, the number of degenerate paths (paths of the same length) increase and so the overall complexity of the aggregate. Secondly, when we build the multilayer transportation network, we add multiple interconnections between each two layers, and we thus increase the degree of nodes that are at the junctions of two layers. It is also important to notice the clear increase of path complexity between the aggregate graph and the different layers taken separately(train, metro, road). The aggregation of layers increases the number of rail degenerated paths from typically one or two to several. The number of degenerated path increases as well but rather slowly with respect to the large number of paths already existing before aggregation.

These effects combined increase the search complexity of a given path in the multilayer transport network and so it increases the difficulty of finding a correct mapping of the sparse trajectories on the graph. These phenomenon also explain why in multimodal transport system using an algorithm that try to find the best match between trajectory(cellular trajectories) and the shape the transport network will usually fails, due to the presence of many degenerate paths.

Framework and Overall Design

Given the multimodal transportation network G and the cellular network Then, in the second phase, (shown in Fig. 5d) the algorithm traverses the skeleton sequence and outputs a sequence of adjacent intersections by completing the sequence (shown in Fig. 5e). Please note that the skeleton sequence searched in the first phase is with equal-length to the given sparse cellular trajectory O, while the intersection sequence outputted in the second phase would be longer than O. Given the frequency of 15 minutes for observations, it is clear that a user would pass through more than one intersection between each 2 consecutive observation points, (e.g. when commuting with metro, it takes around 3 minutes to move between each 2 stations).

Skeleton Sequence Search -Given the sparse cellular trajectory o 0 → o 1 , ... → o M , this phase returns the skeleton sequence of the intersections as

v 0 → v 1 , ... → v M .
P r t (v j ) = P (o t |v j ) × max ∀ vi [P r t-1 (v i ) × T r(v i , v j )] (3)
where P (o t |v j ) is the probability of a user connecting to cell tower of o t when he/she is in the intersection v j and T r(v i , v j ) is the transition probability of moving from node v i to node v j . The parent node is also stored using Eq. 4;

P ar(v j ) = arg max ∀ vi [P r t-1 (v i ) × T r(v i , v j )] (4)
At the end, we find * v M = arg max

tersections * v i , * v i+1 as: * v i → v i1 → v i2 ... → v i k → * v i+1 (6) ⇑ Recovered path
as the complete adjacent sequence. Please note that each two consecutive nodes in the newly obtained sub-sequence are adjacent in multiplex G. In the next section, we will introduce the calculation of probabilities used in our framework.

Core Algorithms

In the previous section we described the general algorithm of mapping cellular trajectories over the multimodal transportation network. The two main probability distributions used in the mapping algorithm, are the HMM transition and emission scores that are estimated in an unsupervised way. This section explains in detail how the two scores are defined and estimated.

Transition Probability

The transition probability T r(v i , v j ) in our mapping algorithm specifies the probability of an individual's moving from hidden state v i at time t-1 to hidden state v j at time t. The transition probability is inferred from the underlying network, i.e. the multilayer transportation network in which each transportation layer has its specific characteristics and properties. Table 1 shows some graph topological properties such as the node degree distribution and the physical edge length distribution in different layers of the multimodal transportation network.

It is crucial to notice that relying on the topological properties of network layers without considering their differences, leads to a biased mapping algorithm in which the observations tend to be mapped over a specific transportation layer.

In addition, taking into account the sparseness of cellular observations, it is a key to authorize transitions between nonadjacent intersections. We propose a transition probability of moving from intersection v i to the intersection v j that is a function of 2 given factors:

1) Edge type and average speed over each edge: each physical edge in the multilayer graph G belongs to a layer. Moreover, only the road layer contains different types of edges (such as highway, principal, local, etc.). We define matrix W where each element of the matrix represents a weight between two nodes if there exists an interconnection between them. The weight of each link is defined as the inverse of average speed that one could have over the corresponding edge.

Table 2 represents the weight according to average speed over the edges of graph G.

W ij =      w ij if v i , v j are adjacent in G 0 otherwise. (7)
2) Edge length: involving edge length in the transition probability, indirectly considers higher probability for the nodes close to each than that for farther ones.

The transition probability between two intersections v i and v j is defined as the inverse of the shortest path cost between v i and v j :

w ij Condition 1/80 Ψ(v i ) = Ψ(v j ) = metro 1/90 Ψ(v i ) = Ψ(v j ) = road (highway) 1/60 Ψ(v i ) = Ψ(v j ) = road (principale) 1/40 Ψ(v i ) = Ψ(v j ) = road (regional) 1/30 Ψ(v i ) = Ψ(v j ) = road (local)
T r(v i , v j ) =   ∑ ∀ (mn) ∈ SP vivj w mn × d(v m , v n )   -1 (8)
where (mn) is an edge between v m and v n belonging to SP vivj , the shortest path between two nodes v i and v j in graph G. The shortest path cost of SP vivj is the sum of distances over each edge (mn) belonging to SP vivj , weighted by

w mn . d(v m , v n )
is the euclidean distance between each two nodes v m and v n .

In earlier studies, the transition probability was quantified based on topological properties of the underlying network which was mainly a road graph.

In [START_REF] Newson | Hidden markov map matching through noise and sparseness[END_REF][START_REF] Thiagarajan | Vtrack: Accurate, energy-aware road traffic delay estimation using mobile phones[END_REF], the transportation network was represented as road segments and transitions were assumed to occur between adjacent road segments. The authors in [START_REF] Thiagarajan | Vtrack: Accurate, energy-aware road traffic delay estimation using mobile phones[END_REF][START_REF] Hummel | Map matching for vehicle guidances[END_REF] considered equal transition probabilities between nodes in the same road segment or nodes between road segments which are adjacent with an intersection. The transition probability in [START_REF] Thiagarajan | Accurate, low-energy trajectory mapping for mobile devices[END_REF] is defined based on the Manhattan distance between the grid cells of the road network. The objective of our proposed transition probability model is to minimize the bias of the mapping algorithm for layers with different topological properties.

Emission Probability

In HMM, at each time step t, there exists an observation o t which in our study is characterized as c t =< lon, lat, r max t > t . The emission score reflects the notion that it is more likely that a particular observation point is observed from a nearby intersection than an intersection farther away [START_REF] Thiagarajan | Vtrack: Accurate, energy-aware road traffic delay estimation using mobile phones[END_REF]. For studies in which GPS data were used as observations [START_REF] Thiagarajan | Vtrack: Accurate, energy-aware road traffic delay estimation using mobile phones[END_REF][START_REF] Newson | Hidden markov map matching through noise and sparseness[END_REF][START_REF] Goh | Online map-matching based on hidden markov model for real-time traffic sensing applications[END_REF], the emission probability score is modeled by a normal distribution that is a function of the euclidean distance between the observation point and the hidden state, and with a standard deviation estimated from sensor errors.

In this work, cellular antenna locations serve as observations; since there is no labeled data available to estimate cellular sensor errors, we build the Voronoi tessellation of cellular antennas in the area of study. In the voronoi network of cellular antennas, each cellular antenna C i is characterized by radius r i which is the maximum distance of the cellular antenna from the corresponding voronoi cell vertices. Our emission score is defined as a decreasing function of the distance between the antenna location and the hidden node (intersection):

P r(o t |v j ) ∝              1.0 if : d tj ≤ r max t ( r max t d tj ) β if : r max t ≤ d tj ≤ τ 0 otherwise. ( 9
)
where d tj = d(o t , v j ) is the euclidean distance between o t and intersection v j , and τ is a threshold corresponding to the maximum distance that a cell phone can be hit by a cellular antenna. τ enforces the constraint that only intersections in the radius of τ from the cellular antenna could be considered as candidate states (nodes).

Evaluation

Dataset for Evaluation

In order to evaluate the proposed algorithm, GPS data are used as ground truth. We collected the cellular trajectories of 10 volunteer participants during one month (Aug-Sept 2014) with their corresponding GPS data. The GPS data were collected with the help of the application "Moves" [START_REF]Moves[END_REF] which was installed on participants' smartphones. The data captured were the sampled phone positions during its movements as well as its activities classified in four different categories: 'Walking', 'Running', 'Cycling' and 'Transport'. Based on this dataset, a set of prepossessing steps were performed in order to extract trajectories mapped over the transport networks.

Furthermore, trajectories whose lengths are shorter than 5 kilometers were filtered out from the database. Given the low sampling rate of cellular data (a data point every 15 minutes), it is not realistic to seek recovering a movement that lasts less than this threshold. The effect of that filter on the dataset distribution could be observed in Fig. 6a and Fig. 6b.

The spatial accuracy needed in order to distinguish a real mobility from noise depends on the distance between two base stations. In order to filter out irrelevant movements, we filtered out all the trajectories under the threshold x th such that P r (X < x th ) = q were P r (X) is the distribution of distance between neighboring antennas, for q = 0.97 as Fig. 7 shows all the neighboring distance is less than 5 kilometers.

As conclusion, we built a dataset of 80 cellular trajectories (sequence of base stations) with their corresponding GPS paths mapped over a multilayer graph G. The multilayer transportation network contains around 16000 nodes and 26000 edges. The users trajectories covered a total distance of 2200 kilometers.

The average number of observation points in each cellular trajectory is 5.55 and the average length of a trajectory is 26.5 kilometers. Fig. 8 shows the coverage area of collected the GPS trajectory dataset.

Evaluation Results and Comparison

Mapping algorithm efficiency

To evaluate our algorithm, the aforementioned labeled dataset was used for test and evaluation. We preformed CT-Mapper to map the cellular trajectories Considering each trajectory as a set of nodes, precision is the fraction of retrieved nodes that belong to the real path. Recall (also known as sensitivity) is the fraction of correct nodes that are retrieved by the algorithm. Moreover, in the evaluation section Root Mean Square Error (RMSE) have been used for two purposes: First, to quantify the overall distance between the obtained result and the ground truth. Second, owing to the considerable spatial noise of cellular observations, RMSE is used to detect matches between 2 points using threshold ϵ .

In this case, if the RMS error between two points is smaller than ϵ , we consider the inferred point as a match. For example, an error of 0.1 kilometers indicates that for each node in the output sequence, the node is considered as a match point if it is within a 0.1 kilometer radius of its corresponding real location. We calculated the four mentioned accuracy results (precision, recall, skeleton and complete sequence similarity score) for a range of fixed allowed RMSE on the obtained mapping results. The similarity scores are the complementary of the Edit distance scores. Fig. 9 represents the result of this evaluation. As Fig. 9 shows, with allowed RMSE of 200 meters, more than 50% of skeleton and complete trajectories can be retrieved. This is remarkable given the sparsity of the coarse grain cellular antenna positions with respect to real user trajectory (average of 5.5 observations per trajectory in the dataset while the average length is 26.5 km). It is important, also, to mention that the frequency of cellular data col-lection is 15 minutes and with higher frequency, higher performance is expected for CT-Mapper. The average similarity score, for a RMSE of 1 kilometer, raises to 80%. In addition, CT-Mapper reaches a recall and a precision of around 80% when a RMSE of 1 kilometer is allowed. In addition to the metrics mentioned above, we compute the Edit distance error not as the number of required edits, but by considering the euclidean distance as the cost of each required edit.

The average of Edit distances for all trajectories in the dataset is 0.79 kilometer.

Comparison with Baseline Algorithms

In this section, the performance of our proposed model is compared with two baseline models. Baseline 1 is a simple model that snaps each observation to the nearest node in the network to find the skeleton and for the second phase, uses least-cost paths between them to retrieve the full path. The result of this baseline model is compare with CT-Mapper in Fig. 10.

To evaluate our transition probability model based on transportation properties as presented in Eq.( 8), we derive Baseline 2, an HMM based baseline model associated with the naive assumption consisting of setting equal probabilities for all outgoing transitions from each node (including self node transition). Under such a model, the transition probability between two nodes v i and v j is represented as:

T r(v i , v j ) =   k i * ∏ n∈Q k n   -1 (10)
where Q = SP vivj -{v i , v j } and k i is the degree of v i . This naive assumption considers all the multilayer network edges on equal footing irrespective of their layer transportation properties .

Using this transition probability model, we build an HMM in the same way as CT-Mapper was developed. We use this model as a baseline algorithm and run it on the test dataset to compare the results with CT-Mapper. We calculate all four performance measures for the baseline models. Fig. 10 compares the performances of the two models with CT-Mapper. As the figures show, there

Multimodality analysis

In the next step of assessing our mapping algorithm, we investigate the accuracy of the mapping algorithm in transportation layer detection. As mentioned in Sec. 8, the complexity of multimodal mapping significantly increases owing to the considerable topological differences between transportation layers.

This issue is dealt with in the proposed transition probability model that seeks minimizing the bias in the mapping algorithm.

We calculate the recall and precision for correct layer detection once for each Fig. 12 shows these measures compared with the baseline algorithm. We have to notice that since each assumption considers specific aspect of network's topological properties, they might introduce a bias in the mapping problem. As it is shown in Fig. 12, the overall recall and precision of correct layer detection is improved in CT-Mapper compared to the baseline algorithm.

Discussion & Conclusion

In this study, we proposed an unsupervised mapping algorithm (CT-Mapper)

to map sparse cellular trajectories over a multimodal transportation network.

We modeled and built the multilayer transportation network of subway, train To validate our transition probability model that better accommodates the complexity of the multimodal transportation network, we compared it with a baseline algorithm that does not take into account the transportation properties of each layer and the results show up to 20% of accuracy improvement of the first over the second. We expect that using a dynamic weight matrix which is compatible with the traffic model at different times of the day, is likely to enhance the mapping results. This issue will be investigated in future studies.

The improvement of accuracy measures of our mapping algorithm by minimizing bias mainly emanating from the multimodality of the transportation network is of great importance which shall be discussed in future contributions. Investigating the possibility of using the proposed mapping algorithm at near real-time (NRT) for traffic monitoring is another direction of further contributions.

Ethics requirements and legal requirements followed during the data collection

Before starting the experiment of collecting cellular data we submitted the experiment protocol to the university ethics committee. Once the experiment started each volunteer signed a legal agreement stipulating that each of them requested access to their cellular localization data for one month (with a sampling interval of 15 min). This request was bond with a legal agreement given us the right to use their data for research purpose only. After one month of retention period after the end of the experiment, the cellular data were directly provided to the volunteers by the telecom operator. They forwarded us their cellular data afterward as well as their GPS traces.

Figure 1 :

 1 Figure 1: A user's trip from Airport CDG to city center of Paris: The road trajectory consists of the sequence of roads that the user passes-by; The GPS trajectory is sampled in minute based frequency; The Cellular trajectory (Full) records each cell tower the user passes-by; The CDR trajectory reports the location of the user's each call during the trip; The Sparse Cellular trajectory is sampled every 15 minutes

•

  We propose an unsupervised trajectory mapping algorithm, namely CT-Mapper, which maps cellular location data over the multimodal transportation network. The multimodal transportation network database was built using different references of geospatial resources. The mapping algorithm is modeled by an HMM where the observations correspond to user cellular trajectories and the hidden states are associated with nodes of the multilayer graph. Transition probability and emission score were modeled based on topological properties of the transportation network and the spatial distribution of antenna base stations. Viterbi decoding algorithm helps reduce the complexity of finding the best match which might enable us to deploy our unsupervised mapping algorithm on large scale mobility data sets in order to estimate multimodal traffic in metropolitan areas.

Figure 2 :

 2 Figure 2: Multilayer representation of different transportation networks

-

  Such a graph is represented as G = (V, E, L, Ψ) where V , E represent the vertices and the edges, L is the set of possible layers. In our study we focused on 3 layers: road, train and subway.

Figure 3 :

 3 Figure 3: Voronoi tessellation of cellular antennas in Ile-de-France

Fig. 2 Definition 2 .Definition 3 .

 223 illustrates how different transportation layers have been aggregated to build a multimodal transportation network. Cellular Network -In this work, we characterize a cellular network as a set of cell towers C = {c 0 , c 1 , ...c P }, where each cell tower c p =< lat, lon, r max > p is characterized by its latitude and longitude in the geographical coordinate system and by r max which is the maximum radius of the voronoi cell the cell belongs to in the voronoi graph built from set C. Please note that the location of each cell tower does not coincide with the location of any intersection in the transportation network i.e., ∀v i ∈ V , ∀c p ∈ C, we have < lat, lon > p ̸ =< lat, lon > i . Sparse Cellular Trajectory -Further we define a sparse cellular trajectory of a user as a sequence of time-stamped locations O = o 0 → o 1 ... → o M , where each time-stamped location o t =< c(t) > refers to the cell tower at time-stamp t the user is observed at. Trajectory Mapping Problem -Given a transportation network G, cell tower network C, and a user sparse cellular trajectory O, our search problem is to find a sequence of intersections v 0 → v 1 ... → v q which the user actually passes by on the transportation network. 3.2. Data Collection and Datasets Three types of data are used in this study: multimodal transportation network data, sparse cellular trajectory data, and GPS trajectory data. The multimodal transportation network data are used to build the multilayer network graph and the mobility model for the mapping algorithms. Cellular trajectories are used for testing while GPS trajectories are used as ground truth and not for training HMM parameters. Sparse Cellular Trajectory Data -In this work we use a new type of cellular trajectory named Sparse Cellular Trajectory. A set of techniques for data collection are used to capture GPRS Tunneling Protocol (GTP) messages from the Cellular Data Network. Packet inspection of GTP-C (GTP control plane) enables us to capture users' localization information at higher frequency than the traditional CDR. The GTP is the tunneling protocol used to carry data traffic over the mobile network (from 2G to LTE) to internet. When a smartphone enables its internet connection (e.g. when it is turned on), a message is sent over the network asking for access. This message contains among other things the identity of the phone and the cell id covering the user. Once the session is established, update messages are sent carrying information like the bearer or the cell id. These messages are triggered when the user moves from a BTS to another or by resource allocation. More detail about the localization

Figure 4 :

 4 Figure 4: Graph Entropy: (A) absolute value of the average entropy of the graph where Savg is the entropy of the real graph and S R is the entropy of the random graph with similar characteristics, (B) is the relative of the average graphs entropy of the paths in the subgraph of the metro, train, road

C

  , we define an algorithm that outputs the most likely path or sequence of intersections given the sequence associated with a user sparse cellular trajectory O. In order to infer the accurate sequence of intersections from the given sparse cellular trajectory, we propose a two-phase unsupervised mapping algorithm: in the first phase, the algorithm searches a sequence of intersections, namely the skeleton sequence, where each two consecutive intersections are not necessarily adjacent (shown in Fig.5c). For this objective we developed an unsupervised Hidden Markov Model inference algorithm that accommodates the sparsity of observations (15 minutes). The hidden states in the HMM are the multimodal graph nodes corresponding to road intersections or metro/train stations. The transition probability in our model takes care of sparsity of observations by permitting transitions between nonadjacent nodes as explained in Sec.4.A . For each observation, a set of hidden states are selected as the candidate states in order to minimize the complexity of search in the graph. Given a sequence of sparse cellular observations, our HMM model outputs the most likely sequence over the multiplex network (there are cases of only 3 or 4 observation points).

  ∀vM P r t (v M ) Then a backtracking iteration using Eq. 5 * v b-1 = P ar( * v b ) for b = [M, ..., 2, 1] (5) retrieves the most likely intersection sequence * v 0 → * v 1 , ... → * v M which produces the most likely path for the sparse cellular trajectory o 0 → o 1 , ... → o M . Sequence * v 0 → * v 1 , ... → * v M serves as input for the next phase to retrieve the adjacent sequence of intersections for the given sparse cellular trajectory.

Figure 5 :

 5 Figure 5: An illustration of different phases of mapping algorithm. The Blue line in the Fig. 5(a) is the real GPS trajectory of a user and given a sequence of 5 antenna base stations with the frequency of 15 min, the mapping algorithm can retrieve the pink line in Fig.5(e)

Figure 6 :

 6 Figure 6: Time distribution and distance distribution

Figure 7 :

 7 Figure 7: Neighboring cell distance distribution

Figure 8 :

 8 Figure 8: The coverage area of GPS data collected is shown in yellow on the map of Paris and region

Figure 9 :

 9 Figure 9: Result evaluation

Figure 10 :

 10 Figure 10: Up-left: Precision, up-right: Recall, bottom-left is Edit-based similarity scores and bottom-left is the skeleton similarity score

Figure 11 :Figure 12 :

 1112 Figure 11: Sequence Edist Distance

  and road layers for the Ile-de-France metropolitan area. The multilayer transportation network contains around 16000 nodes and 26000 edges. Investigating the complexity of the multilayer transportation graph, a transition probability model leveraging the transportation layer type and topological properties was estimated and used in an unsupervised HMM-based mapping algorithm. We carried experiments on a test dataset of 80 real multimodal trajectories collected from 10 participants during one month (Aug-Sept 2014) to evaluate our algorithm. Considering the sparsity of cellular observations (with a frequency of 15 minutes), the percentage of retrieved paths of smartphone users is notable.

Table 1 :

 1 Mapping Problem in the Collected Datasets Different transportation networks with their properties

		Number of		Avg.
		Node Edge	Degree Edge	Reference
					length
	Subway	303	356	2.35	0.757	OSM
	Train	241	244	2.025	3.07	OSM
	Tramway 146	140	1.918	0.71	OSM
	Road	14798 22276 3.01	1.34	IGN

Table 2 :

 2 Edge classification and weights for multilayer transportation network G.

In the rest of paper, we use the term "cellular trajectory" and "sparse cellular trajectory" interchangeably.

Acknowledgment

This research received funding from the Pierre and Marie Curie University through their PhD founding program.


Francesco Bonaldi

Giuseppe Geymonat

Françoise Krasucki

Marina Vidrascu

Mathematical and numerical modeling of plate dynamics with rotational inertia

Keywords: Plates, Kirchhoff-Love, Rotational inertia, FreeFEM, HCT, Newmark. 2010 Mathematics Subject Classification. 74K20, 65M60, 35Q74

We give a presentation of the mathematical and numerical treatment of plate dynamics problems including rotational inertia. The presence of rotational inertia in the equation of motion makes the study of such problems interesting. We employ HCT finite elements for space discretization and the Newmark method for time discretization in FreeFEM++, and test such methods in some significant cases: a circular plate clamped all over its lateral surface, a rectangular plate simply supported all over its lateral surface, and an L-shaped clamped plate.

Introduction

In [START_REF] Bonaldi | An asymptotic plate model for magneto-electro-thermo-elastic sensors and actuators[END_REF] the equations of plate models for magneto-electro-thermo-elastic sensors and actuators have been deduced by an asymptotic development with respect to the thickness 2e ą 0 of the plate. A peculiar feature of the different models concerns the flexural problem, governing the evolution of the transversal displacement w of the plate, which occupies the closure of the three-dimensional domain Ω ˆp´e, eq in its reference configuration and whose kinematics is of Kirchhoff-Love type. This problem is uncoupled from the membrane problem, it takes into account an inertia effect involving the mean curvature of the deformed middle surface, referred to as rotational inertia, and the only influence of magneto-electro-thermo-elastic behavior of the material appears in the coefficients A " pA αβστ q of the (symmetric) moment tensor Mptq " pM αβ ptqq " ´2e 3 3 A∇∇wptq, p∇∇wptqq αβ " B αβ wptq (let us point out that the fourthorder tensor A " pA αβστ q is symmetric and strongly elliptic). The transverse displacement w of the plate is solution of the following evolution equation:

2eρ : w ´2e 3 3 ρ : w ´div div M " f `div m in Ω ˆp0, T q, (1)
equipped with initial conditions wp0q " w 0 , 9 wp0q " w 1 in Ω

and with suitable boundary conditions, where ρ ą 0 is the mass density of the plate, supposed to be homogeneous. If one supposes, for the sake of simplicity, that the three-dimensional body is subjected to time-dependent volume forces p0, 0, Fptqq and surface loads pG 1 ptq, G 2 ptq, 0q on the upper and lower faces Ω ˆt˘eu, then the source terms in the right-hand side of (1) have the following expressions f ptq "

ż `e ´e Fptq dx 3 , m α ptq " epG ὰ ptq ´Gά ptqq, α P t1, 2u.

The term accounting for rotational inertia in (1) is ´2e 3 3 ρ : w, and it derives from the expression of the volume kinetic energy kptq for Kirchhoff-Love plates subjected to pure bending:

kptq " 1 2 ˆ2eρ 9 w 2 `2e 3 3 ρ|∇ 9 w| 2 ˙.
The reason for the word rotational is that the quantity |∇w| gives a measure of the rotation of the material fibers orthogonal to the middle surface of the plate, which remain straight and perpendicular to this surface after deformation in Kirchhoff-Love kinematics.

The influence of rotational inertia on the lateral vibrations of linearly elastic bars was considered by Lord Rayleigh [START_REF] Baron Rayleigh | The Theory of Sound[END_REF], sect. 186; the extension of this analysis to the flexural motion of isotropic elastic plates was carried out for the first time in 1951 in a seminal paper by R. D. Mindlin [START_REF] Mindlin | Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates[END_REF] where also the effect of transverse shear deformation is taken into account. Later on, an evolution model for plates with rotational inertia was deduced by A. Raoult [START_REF] Raoult | Construction d'un modèle d'évolution de plaques avec terme d'inertie de rotation[END_REF] using the asymptotic expansion method. The influence of such effects has also been considered in the case of flexural motion of large amplitude. Analogously to the case of beams, where it can be worth to account for rotational inertia of the cross section in case the slenderness of the beam (i.e., the ratio of the cross section diameter to the length of the beam) is lower than one, but not very close to zero, in our case rotational inertia of the material fibers orthogonal to the middle surface may be important if the plate thinness (i.e., the ratio of the thickness to the diameter of the plate) is lower than one, but not very close to zero. In the following, to point out the influence of the thickness parameter, we carry out our analysis on equation (1) divided by 2e.

Our goal, after providing a result of well-posedness for the problem, is to test a numerical method involving C 1 elements for what concerns space discretization and Newmark's midpoint method concerning time discretization. The presence of the rotational inertia requires an appropriate choice of the functional spaces, presented in sect. 1, where we sketch the proof of the well-posedness using the Faedo-Galerkin method. In sect. 2 we perform, in the appropriate Sobolev spaces, the error analysis of a finite element spatial discretization and of a Newmark-type discrete time approximation. Our choice allows the application of the methods developed for linear second-order evolution equations. Based on a continuoustime Galerkin method (see e.g. [START_REF] Baker | Error estimates for finite element methods for second order hyperbolic equations[END_REF] or [START_REF] Dupont | L 2 -Estimates for Galerkin Methods for Second Order Hyperbolic Equations[END_REF]) we can infer optimal error estimates, and then couple with the error estimates for the time discretization of Newmark type (see e.g. [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF]). Note that a conforming space discretization for plates requires C 1 elements; we choose HCT elements (see e.g. [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]). This theoretical numerical analysis is complemented with numerical tests performed with FreeFEM++.

Let us remark that the choice of C 1 elements is not very common in practice, as these elements are rather expensive. Mixed and hybrid plate elements were extensively studied. Among the elements widely used in the engineering community we can mention the DKT triangle [START_REF] Batoz | A study of three-node triangular plate bending elements[END_REF]. For a general review of plate elements see for example [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF], [START_REF] Hughes | The finite element method: linear static and dynamic finite element analysis[END_REF] and references therein. We can also mention the more modern hybrid high-order methods [START_REF] Di Pietro | An Arbitrary-Order and Compact-Stencil Discretization of Diffusion on General Meshes Based on Local Reconstruction Operators[END_REF][START_REF] Di Pietro | A hybrid high-order locking-free method for linear elasticity on general meshes[END_REF]: they will be studied in a forthcoming paper.

Existence and uniqueness

We start by rewriting (1) as

ρ : w ´e2 3 ρ : w ´e2 3 div div A∇∇w " 1 2e pf `div mq. (4)
The weak formulation of the problem depends on boundary conditions. For this, we first define the pivot space

H " tu P H 1 pΩq : u " 0 on Γ 0 u (5)
and the bilinear form b : H ˆH Ñ R given by bpu, vq "

ż Ω ˆρuv `e2 3 ρ∇u ¨∇v ˙dΩ, @u, v P H. (6)
Notice that bp¨, ¨q defines a scalar product in H whose associated norm | ¨|b is equivalent to the usual Sobolev norm. Let V be a Hilbert-Sobolev space such that

piq V Ď H 2 pΩq X H, piiq the embedding V ãÑ H is compact.
Thanks to piq, we can endow V with the Hessian L 2 -norm given by }v} "

ˆżΩ ∇∇v : ∇∇v dΩ ˙1{2 . (7)
Let us define the bilinear form a : V ˆV Ñ R by apu, vq " e 2 3

ż

Ω A∇∇u : ∇∇v dΩ, @u, v P V. (8)
By the symmetry and ellipticity properties of A, there exist two positive constants A ´and A `such that, for all x P Ω,

A ´|U| 2 ď AU :U ď A `|U| 2 , @U P Symp2q.
Hence, ap¨, ¨q is symmetric and V-elliptic: apv, vq ě e 2 3 A ´}v} 2 for any v P V . In the numerical examples, we consider essentially the following situations:

pBCq 1 : $ ' & ' % 2e 3 3 ρ B n :
w `div M ¨n `Bτ pMn ¨τ q " ´m ¨n on Γ 1 ˆp0, T q, Mn ¨n " 0 on Γ 1 ˆp0, T q, w " 0, B n w " 0 on Γ 0 ˆp0, T q, pBCq 2 :

$ ' & ' % 2e 3
3 ρ B n : w `div M ¨n `Bτ pMn ¨τ q " ´m ¨n on Γ 1 ˆp0, T q, Mn ¨n " 0 on BΩ ˆp0, T q, w " 0 on Γ 0 ˆp0, T q.

Boundary conditions pBCq 1 refer to a plate clamped on Γ 0 ; in this case we choose

V " V 1 " tu P H 2 pΩq : u " B n u " 0 on Γ 0 u. (9)
Boundary conditions pBCq 2 feature a plate simply supported on Γ 0 , in which case we choose

V " V 2 " tu P H 2 pΩq : u " 0 on Γ 0 u " H 2 pΩq X H. (10)
In both cases, under the general hypothesis of anisotropic linearly elastic behavior, the weak formulation of the problem has the following aspect:

$ ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' %
For any fixed t P p0, T q, find wptq P V such that

ż Ω ˆρ : wptqv `e2 3 ρ∇ : wptq ¨∇v `e2 3
A∇∇wptq : ∇∇v ˙dΩ " 1 2e
ż Ω pf ptqv ´mptq ¨∇vq dΩ, for all v P V, with initial conditions wp0q " w 0 , 9 wp0q " w 1 ,

where

V " V 1 or V " V 2 .
Remark 1.1. Notice that, by taking volume and surface loads constant in space in (3), we have f ptq " 2eFptq and m α ptq " epG ὰ ptq´G ά ptqq, so that the right-hand sides of (4) and of (11) are independent of the thickness.

In order to show that problem [START_REF] Dupont | L 2 -Estimates for Galerkin Methods for Second Order Hyperbolic Equations[END_REF] is well-posed, we identify the time-dependent linear form on H L t pvq " 1 2e

ż Ω pf ptqv ´mptq ¨∇vq dΩ, with the scalar product of an element F ptq P H (for 0 ă t ă T ) with v P H. This can be accomplished via the following problem:

Find F ptq P H such that bpF ptq, vq " L t pvq, @v P H.

Provided f ptq P L 2 pΩq and mptq P L 2 pΩq for 0 ă t ă T , problem ( 12) is well-posed by the Lax-Milgram Lemma. As for problem [START_REF] Dupont | L 2 -Estimates for Galerkin Methods for Second Order Hyperbolic Equations[END_REF], we assume f P L 2 p0, T ; L 2 pΩqq and m P L 2 p0, T ; L 2 pΩqq, so that F P L 2 p0, T ; Hq; the formulation of the problem reads then:

Given w 0 P V , w 1 P H and F P L 2 p0, T ; Hq, find a function w such that w P C 0 pr0, T s; V q X C 1 pr0, T s; Hq and, for all v P V, d 2 dt 2 bpwptq, vq `apwptq, vq " bpF ptq, vq, wp0q " w 0 , dw dt p0q " w 1 .

(

) 13
Theorem 1.2. Let T ą 0 be fixed, w 0 P V , w 1 P H and F P L 2 p0, T ; Hq.

(i) There exists a unique function w P C 0 pr0, T s; V q X C 1 pr0, T s; Hq satisfying (13).

(ii) For all t P r0, T s, the function w satisfies the following energy equation:

| 9 wptq| 2 b `apwptq, wptqq " |w 1 | 2 b `apwp0q, wp0qq `2 ż t 0 bpF psq, 9 wpsqq ds.
(iii) There exists a constant c " cpΩ, Aq such that

}w} C 0 pr0,T s;V q `}w} C 1 pr0,T s;Hq ď ? T e cpΩ, Aq ´}w 0 }`|w 1 | b `}F } L 2 p0,T ;Hq ¯.
For a proof of statements (i) and (iii) of Theorem 1.2, of which we present hereinafter some essential parts, see e.g. [START_REF] Allaire | Numerical Analysis and Optimization[END_REF] and [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF] (see also [START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology[END_REF], for a general treatment of evolution equations); the proof of statement (ii) can be found in [START_REF] Lions | Non-Homogeneous Boundary Value Problems and Applications[END_REF]Chap. 3,Sect. 8] where a very general situation is considered. The uniqueness of the solution is proven in a standard way. With the choice ( 9) or (10) of the space V , and (5) for the space H, endowed with scalar products ( 8) and ( 6) respectively, the Faedo-Galerkin approximation method can be used. The compactness of the embedding V ãÑ H implies that there exists an increasing sequence of eigenvalues 0 ă λ 1 ď λ 2 ď . . . ď λ i ď . . . and a Hilbert basis tg i u, orthonormal in H and orthogonal in V , of associated eigenvectors verifying @v P V, apg i , vq " λ i bpg i , vq (see e.g. [?]). As for existence, the subspace V m of V generated by the first m eigenvectors g 1 , . . . , g m is introduced; let w m : t P r0, T s Þ Ñ w m ptq P V m be solu-tion to the following (well-posed) second-order system of differential equations:

@v m P V m , d 2 dt 2 bpw m ptq, vq `apw m ptq, vq " bpF ptq, vq, (14)
w m p0q " w 0,m " m ÿ i"1 bpw 0 , g i qg i , dw m dt p0q " w 1,m " m ÿ i"1 bpw 1 , g i qg i . ( 15
)
One can show that tw m u mPN is a Cauchy sequence in spaces C 0 pr0, T s; V q and C 1 pr0, T s; Hq and thus it converges to a function w in such spaces. Finally, by a density argument, it is proven that w is solution to [START_REF] Hecht | New development in FreeFem++[END_REF].

2 Numerical analysis

Semi-discrete problem

Let V h Ă V denote a subspace of V of dimension I " Iphq and consider the following semi-discrete problem: given w 0,h P V h and w 1,h P V h , find the solution w h : t P r0, T s Þ Ñ w h ptq P V h to the following system of ordinary differential equations:

@v h P V h , d 2 dt 2 bpw h ptq, v h q `apw h ptq, v h q " bpF ptq, v h q, w h p0q " w 0,h , dw h dt p0q " w 1,h . (16)
We now introduce a basis tϕ i u 1ďiďI of V h and denote the time-dependent components of w h in this basis by ξ j ptq, j " 1, . . . , I. Analogously, we denote the components in the same basis of w 0,h and w 1,h respectively by ξ 0,j and ξ 1,j . Finally, we set χ hj ptq " bpF ptq, ϕ j q. Then ( 16) reads

# M h : ξptq `Kh ξptq " χ h ptq, ξp0q " ξ 0 , 9 ξp0q " ξ 1 , (17)
with self-explanatory notation. Matrices M h and K h are respectively the mass matrix and the stiffness matrix. Their coefficients are pM h q ij " bpϕ i , ϕ j q and pK h q ij " apϕ i , ϕ j q, 1 ď i, j ď I.

In the following, we drop the index h when there is no ambiguity.

Remark 2.1. The choice of the scalar product ( 6) in H ensures that M is positive definite.

Problem [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF] is the semi-discrete counterpart of (13). Under certain in-time regularity hypotheses on w and adapting the arguments of [START_REF] Baker | Error estimates for finite element methods for second order hyperbolic equations[END_REF], [START_REF] Dupont | L 2 -Estimates for Galerkin Methods for Second Order Hyperbolic Equations[END_REF], it is possible to give an estimate of the error w h ptq ´wptq. For this, we introduce the elliptic projection operator Π h , a linear and continuous operator mapping u P V onto Π h u P V h defined by @v h P V h , apΠ h u ´u, v h q " 0.

We have the following results.

Theorem 2.2. Let T ą 0 be fixed and assume that the solution w of (13) verifies w P C 2 pr0, T s; V q. Then there exists a constant C " CpAq independent of h such that, for any t P r0, T s,

}w h ptq ´wptq} `ˇˇˇd w h dt ptq ´dw dt ptq ˇˇˇb ď 1 e CpAq ˆ}w 0,h ´Πh w 0 } `|w 1,h ´Πh w 1 | b `}pI ´Πh qwptq} `ˇˇˇp I ´Πh q dw dt ptq ˇˇˇb `ż t 0 ˇˇˇp I ´Πh q d 2 w dt 2 psq ˇˇˇb ds ˙.
Theorem 2.3. Under the assumptions of Theorem 2.2 let the following approximation hypotheses be satisfied: w h " w in C 0 pr0, T s; V q and in C 1 pr0, T s; Hq.

@v P V, lim hÑ0 inf v h PV h }v ´vh } " 0, (19)
lim hÑ0 }w 0,h ´w0 } " 0, (20)

Time discretization: Newmark's midpoint method

A typical time discretization method widely used for differential systems of the form ( 17) is the Newmark method (see, e.g., [START_REF] Géradin | Mechanical Vibrations: Theory and Application to Structural Dynamics[END_REF]). In particular, we will use a midpoint approximation. Let N P N; we introduce a time-step ∆t " T N and a uniform mesh of the interval r0, T s defined by nodes t n " n∆t, 0 ď n ď N . We want to compute an approximation ´ξn , 9

ξ n ¯of the pair ´ξpt n q, 9 ξpt n q ¯. The discretization of (17) based on Newmark's midpoint method can be rewritten as

$ ' & ' % ˆM `∆t 2 4 K ˙ξn`1 2 " Mξ n `∆t 2 M 9 ξ n `∆t 2 4 χ n`1 2 , 0 ď n ď N ´1, ξ 0 " ξ 0 , 9 ξ 0 " ξ 1 (22) with ξ n`1 2 " 1 2
`ξn`1 `ξn ˘and χ n`1 2 " 1 2 pχpt n`1 q `χpt n qq. Upon solving (22) for ξ n`1 2 , the update rules are

9 ξ n`1 " 4 ∆t ´ξn`1 2 ´ξn ¯´9 ξ n , ξ n`1 " 2ξ n`1 2 ´ξn .
In this situation, Raviart and Thomas [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF] have proved the following error estimate.

Theorem 2.4. Let T ą 0 be fixed. Then the solution tw n h P V h , 0 ď n ď N u given by (22) verifies the following estimates:

(i) if w P C 2 pr0, T s; V q X C 3 pr0, T s; Hq, |w n h ´wpt n q| b ď C ˆ|w 0,h ´Πh w 0 | b `|w 1,h ´Πh w 1 | b |pI ´Πh qwpt n q| b `ż tn 0 ˆˇˇˇp I ´Πh q d 2 w dt 2 psq ˇˇˇb `∆t ˇˇˇd 3 w dt 3 psq ˇˇˇb ˙ds ˙,
where C " CpT q is independent of h, ∆t and w;

(ii) if w P C 2 pr0, T s; V q X C 4 pr0, T s; Hq,

|w n h ´wpt n q| b ď C ˆ|w 0,h ´Πh w 0 | b `|w 1,h ´Πh w 1 | b |pI ´Πh qwpt n q| b `ż tn 0 ˆˇˇˇp I ´Πh q d 2 w dt 2 psq ˇˇˇb `∆t 2 ˇˇˇd 4 w dt 4 psq ˇˇˇb ˙ds ˙,
where C " CpT q is independent of h, ∆t and w.

Space discretization: HCT elements

The solution of ( 13) belongs to C 0 pr0, T s; V q X C 1 pr0, T s; Hq, thus we are led to select finite elements of class C 1 ; in particular, we will use HCT elements.

Let T h denote a regular mesh in the sense of Ciarlet [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF][START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF] of the domain Ω, K P T h the typical element of T h and X h a finite element space. Moreover, let P K " tv h|K : v h P X h u. We recall the following result of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF].

Theorem 2.5. Assume that the inclusions P K Ă H 2 pKq for all K P T h and X h Ă C 1 pΩq hold. Then the following inclusions hold:

X h Ă H 2 pΩq, X oh " tv h P X h : v h " 0 on Γ 0 u Ă V 2 , X ooh " tv h P X h : v h " B n v h " 0 on Γ 0 u Ă V 1 .
Remark 2.6. The choices V h " X ooh for pBCq 1 and V h " X oh for pBCq 2 ensure hypotheses [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF] to (21) of Theorem 2.3 to be satisfied [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF].

Finite elements of class C 1 are rather complicated and time-consuming, and are not used too often in practical applications. We choose to start our experiments with such elements because the theory described is valid for conforming approximations. In this context, the HCT is one of the simplest C 1 elements (Fig. 1(a)). The set of degrees of freedom (twelve in total) is given by the values of a function, as well as of its partial derivatives, at the three vertices and by the values of its normal derivatives at the midpoints of the sides. From an internal viewpoint, the HCT element is a composite element: a typical triangle K P T h is split into three sub-triangles K i (i " 1, 2, 3), the internal node usually corresponding to the barycenter of K (Fig. 1(b)). A polynomial of degree three is defined on each sub-triangle, so that the space P K is given by P K " tp P C 1 pKq : p |K i P P 3 pK i q, 1 ď i ď 3u.

The condition p P C 1 pKq is realized by requiring the continuity of the three polynomial expansions and of their gradients at the barycenter (marked by a black square in Fig. 1(b)), and the continuity of their normal derivatives at the midpoints of the internal sides. Thus, the HCT element is a C 1 element as a whole, in the sense that a function and its first derivatives are continuous across the edges of any two adjacent elements of T h .

Let us now provisionally focus on the static counterpart of (13), namely,

@v P V, apw, vq " Lpvq, (23)
where we have neglected time dependence. Notice that in this case the pivot space H is L 2 pΩq. The discrete version of (23) reads

I ÿ j"1
apϕ i , ϕ j qξ j " Lpϕ i q, 1 ď i ď I.

The following theorem [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] yields an estimate of the error between w h and w when HCT elements are employed in space discretization for (23).

Theorem 2.7. If the exact solution w P V of (23) is also in the space H 4 pΩq, then there exists a constant C ą 0 independent of h such that

}w ´wh } H 2 pΩq ď Ch 2 |w| H 4 pΩq . ( 24
)
Remark 2.8 (Implementation issues). In finite element methods, a quadrature scheme is needed to compute the coefficients apϕ i , ϕ j q and Lpϕ i q, thereby resulting in an approximated bilinear form a h p¨, ¨q and in an approximated linear form L h p¨q. Generally, integration over a mesh element is performed using a quadrature scheme for which all nodes are situated at the interior of the element. However, in the HCT case, a mesh element features internal interfaces between any two sub-triangles; at these interfaces, the continuity of second partial derivatives is not guaranteed. Hence, one should use a quadrature scheme that avoids nodes on any such interface. A solution is to integrate on each sub-triangle and then sum up the three contributions. Moreover, in order to apply the first Strang lemma [?, 7], we need that the bilinear form a h be uniformly V h -elliptic, i.e. Dα ą 0 : @v h P V h , a h pv h , v h q ě α}v h } 2 . In our situation, if the space P K contains polynomials of degree at most k, a sufficient condition for a h to satisfy this property is that the quadrature scheme be exact for polynomials of degree 2k ´4 at least. This means that on each K i Ă K, a quadrature formula exact at least for polynomials of degree two has to be used.

Numerical simulations

For what concerns the numerical treatment of problem ( 1)-( 2), we will restrict our attention to an isotropic behavior; hence, the constitutive equation yielding M is

M " ´D pp1 ´νq∇∇u `ν u Iq , with D " 2e 3 E 3p1 ´ν2 q (25)
the flexural rigidity of the plate, E and ν being respectively Young's modulus and Poisson's ratio of the material. We perform numerical tests using the software package FreeFEM++ 3.42 (see [START_REF] Hecht | New development in FreeFem++[END_REF]), in which HCT elements have been implemented along with an adequate quadrature formula for their use. We will consider two situations:

p1q Ω " tpx, yq P R 2 : x 2 `y2 ă R 2 u, Γ 0 " BΩ and pBCq 1 , i.e. a circular plate of radius R clamped all over the lateral surface;

p2q Ω " p0, aq ˆp0, bq with a, b ą 0, Γ 0 " BΩ and pBCq 2 , i.e. a rectangular plate simply supported all over the lateral surface. In both situations we assume m " 0, and in order to get coherent results we fix once and for all the following set of data: R " 5 cm, a " 6 cm, b " 8 cm, e " 1 mm ρ " 5600 kg{m 3 , E " 136 GP a, ν " 0.3, pi.e. D " 99.63 N ¨mq.

(26)

Statics

In each of the following test cases, we consider nested meshes in order to determine the behavior of the relative error in H 2 pΩq-norm with respect to the meshsize h; of course, mesh refinement is uniform.

Boundary conditions pBCq 1

The problem formulation reads in this case

# D w " f 0 in Ω, w " 0, B n w " 0 on BΩ, (27)
with f 0 a constant; the closed-form solution of ( 27) is given in this case by wpx, yq "

f 0 64D pR 2 ´px 2 `y2 qq 2 . ( 28
)
Table 1. Variation of the H 2 -norm of the relative error with the number of finite elements for w given by (28), along with the corresponding number of degrees of freedom.

Relative error in H 2 -norm 0.75% 0.2% 0.06% 0.01% Table 1 shows that the variation of the H 2 -norm of the error with respect to the number of finite elements (and thus to the meshsize) is in agreement with error estimate (24), i.e. quadratic. A different simulation where the exact solution is non-polynomial, namely, wpx, yq " f 0 64D pR 2 ´px 2 `y2 qq 2 sinpaxq, gives similar results.

Boundary conditions pBCq 2

The problem formulation is

$ & % D w " f 0 sin ´π a x ¯sin ´π b y ¯in Ω, w " 0, Mn ¨n " 0 on BΩ, (29)
The closed-form solution is given by wpx, yq " W 0 sin ´π a x ¯sin ´π b y ¯, W 0 "

f 0 π 4 D ˆ1 a 2 `1 b 2 ˙´2 . ( 30
)
The variation of the relative error is shown in Table 2, and once more the convergence is quadratic.

Table 2. Variation of the relative error in H 2 -norm with the number of finite elements for w given by (30), along with the corresponding number of degrees of freedom.

L-shaped clamped plate

To further test the effectiveness of HCT elements, we also performed a numerical simulation in the case of an L-shaped clamped plate, i.e. with boundary conditions pBCq 1 . Table 3 shows the variation of the relative error in H 2 -norm. Since the closed-form solution is not known in this case, comparison has been made with a numerical solution obtained on a fine adapted mesh, consisting of 137232 elements (meshsize 0.08 mm). Computations have been carried out, again, with nested meshes, taking into account the singularity at the inside corner with a suitablyrefined initial mesh. The corvergence rate turns out to be slightly slower than quadratic (1.74). Since the solution has very large variations around the inside corner, adaptive mesh refinement is necessary; however, such numerical aspects are beyond the scope of this paper.

Table 3. Variation of the H2 -norm of the relative error with the number of finite elements for w given by (28), along with the corresponding number of degrees of freedom.

Dynamics

In order to test the accuracy of the Newmark midpoint method combined with HCT elements, we consider the time evolution of the transverse displacement t Þ Ñ w h px 0 , y 0 ; tq, where px 0 , y 0 q is the center of the plate; we consider nested meshes in the two cases pBCq 1 and pBCq 2 . The data set is given in (26). We consider an exact solution of the form wpx, y; tq " kpx, yq sinpβtq, where β is a constant 1 , and we take into account the evolution of the relative error

ε h ptq " $ ' & ' %
ˇˇˇw h px 0 , y 0 ; tq ´wpx 0 , y 0 ; tq wpx 0 , y 0 ; tq ˇˇˇi f t ą 0, wpx 0 , y 0 ; tq ‰ 0,

0 if t " 0. (31)
Since convergence occurs only for 0 ă t ď T where T is fixed, we precise in each example the value of T .

Influence of space discretization

In this section we fix once and for all the time-step to ∆t " 0.01 s, and we point out the influence of space discretization.

Boundary conditions p p pBCq q q 1

The function wpx, y; tq "

f 0 64D pR 2 ´px 2 `y2 qq 2 sinpβtq, (32)
with f 0 a constant, is solution to the problem

$ ' ' ' ' ' & ' ' ' ' ' % 2eρ : w ´2e 3 3 ρ : w `D w " f in Ω ˆp0, T q, w 0 px, yq " 0, w 1 px, yq " f 0 β 64D pR 2 ´px 2 `y2 qq 2 in Ω,
w " 0, B n w " 0 on BΩ, with f px, y; tq " ´β2 f 0 eρ 96D ´3pR 2 ´px 2 `y2 qq 2 `8e 2 pR 2 ´2px 2 `y2 qq ¯sinpβtq `f0 sinpβtq.

The evolution of the error given by (31) is shown in Fig. 3(a), and it reflects the expected behavior: the error evolution is attenuated upon refining the mesh.

Boundary conditions p p pBCq q q 2

The function wpx, y; tq " W 0 sin ´π a x ¯sin ´π b y ¯sinpβtq, The evolution of the relative error given by ( 31) is shown in Fig. 3(b). Let us remark that the behavior of the error corresponding to the finest mesh (represented by a red line) is in this case almost imperceptible, inasmuch as it is very close to zero. Also, notice that the convergence is in this case remarkably faster than in Fig. 3(a); this can be related to the fact that no approximation error concerning the domain geometry is committed, unlike the case of pBCq 1 , where the domain under consideration is a circle.

Influence of time discretization

In this section we consider two different values of the time-step ∆t (namely, 0.01 s and 0.05 s) in order to point out the influence of this parameter.

In the case of pBCq 2 , when the source term vanishes (f " 0), the closedform solution to the dynamic problem can be obtained by separation of variables. Indeed, one has the Fourier development where

ω mn " π 2 ˆm2 a 2 `n2 b 2 ˙g f f e D 2eρ `2 3 π 2 e 3 ρ ´m2 a 2 `n2 b 2 ¯,
and coefficients g 1 mn and g 2 mn are determined by initial conditions. We consider then the following situation: Figure 3. Evolution of the relative error t Þ Ñ ε h ptq with px 0 , y 0 q " p0, 0q for pBCq 1 (a) (continuous line: 27 elements; dashed line: 108 elements; red line: 432 elements), and with px 0 , y 0 q " pa{2, b{2q for pBCq 2 (b) (continuous line: 12 elements; dashed line: 48 elements; red line: 192 elements). In both cases, β " 1 s ´1 and T " 3.5 s. In both cases, all three evolutions become slightly irregular (more remarkably in case (b), with the largest number of elements) when t is close to π; indeed, gpπq " 0 for β " 1 s ´1.

in which case the exact solution is wpx, y; tq " α ω 11 sin ´π a x ¯sin ´π b y ¯sinpω 11 tq.

Given that the exact solution is p2π{ω 11 q-periodic in time, in order to test our numerical method we consider a reasonable value of ω 11 ; say, ω 11 " 10 s ´1. The evolution of the relative error, given by (31), is shown in Fig. 4(a) for a time-step ∆t " 0.05 s and in Fig. 4(b) for a time-step ∆t " 0.01 s. For ∆t " 0.05 s, the obtained behavior is unexpected: mesh refinement results in an amplification of the relative error ; decreasing the time-step to ∆t " 0.01 s yields the expected behavior. Indeed, note that the period τ corresponding to ω 11 " 10 s ´1 is approximately equal to 0.63 s, so that the ratio ∆t{τ is about 0.08 for ∆t " 0.05 s and about 0.01 for ∆t " 0.01 s. The influence of such ratio is well-known in the case of Newmark's midpoint method (see, e.g., [START_REF] Géradin | Mechanical Vibrations: Theory and Application to Structural Dynamics[END_REF]). Figure 4. Evolution of the relative error t Þ Ñ ε h ptq , with px 0 , y 0 q " pa{2, b{2q, for a vanishing initial displacement and a nonzero initial velocity, for a time-step ∆t " 0.05 s (a) and a time-step ∆t " 0.01 s (b). Continuous line: 28 elements; dashed line: 112 elements; red line: 448 elements. In both cases, α " 1 cm{s and T " 0.7 s.

Concluding remarks

We have pointed out the importance of our choice of spaces V and H, required for the treatment of the rotational inertia term, not only in the proof of the problem's well-posedness, but also in error estimates concerning our numerical method. Our simulations, performed using FreeFEM++ 3.42 show that the presence of the rotational inertia term does not affect the efficiency of the Newmark time discretization method combined with conforming finite elements such as HCT elements.

We have tested the implementation of such elements in three cases: clamped circular plate, simply supported rectangular plate and L-shaped clamped plate. In all these cases, our numerical experiments validate the convergence rate of the error predicted by the theory. As is well-known, HCT elements are computationally expensive; it would then be of interest to use nonconforming space discretization methods (mixed or hybrid) [START_REF] Boffi | Mixed finite element methods and applications[END_REF][START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF], such as HHO methods [START_REF] Di Pietro | An Arbitrary-Order and Compact-Stencil Discretization of Diffusion on General Meshes Based on Local Reconstruction Operators[END_REF][START_REF] Di Pietro | A hybrid high-order locking-free method for linear elasticity on general meshes[END_REF]. The application of such methods to plate problems seems particularly interesting and will be carried out in a forthcoming work.

Figure 1 .

 1 Figure 1. The HCT element (a) and the HCT element in detail (b).

Figure 2 .

 2 Figure 2. L-shaped domain.

w

  0 px, yq " 0, and w 1 px, yq " α sin ´π a

F

We obtained similar results for time dependencies of the form t

and arctanpβtq.

Acknowledgments. We thank Frédéric Hecht for his support and collaboration.

This work has been partially supported by the French Agence Nationale de la Recherche (ANR) under grant ARAMIS (Projet «Blanc», N. ANR 12 BS01-0021)