Juan Antonio

Cordero Fuertes

Matthias Philipp

Emmanuel Baccelli

Routing Across Wired and Wireless Mesh Networks: Experimental Compound Internetworking with OSPF

As wireless mesh networks are deployed, a new concept emerges: compound internetworks, i.e., internetworks that contain both wired networks and wireless mesh networks. Routing is one of the key challenges that arises in compound internetworks: indeed, while specific routing protocols are typically used for wired networks on one hand, and for wireless mesh networks on the other hand, it has been observed that operating a single routing protocol to manage a compound internetwork as a whole brings several advantages. In this realm, the IETF has thus standardized protocol extensions to Open Shortest Path First (OSPF, the routing protocol used by more than 50 % of the wired routers in today's Internet), enabling OSPF to operate simultaneously on wired networks, and on wireless mesh or moderately mobile ad hoc networks (MANETs). This paper evaluates the performance of OSPF coupled with such a protocol extension for MANETs on an experimental compound internetwork testbed. This paper reports on experiments carried out with OSPF operating simultaneously over Ethernet and 802.11b. Despite the limitations of the testbed, these experiments provide both a proof-of-concept and complementary results compared to prior work in the domain, which was mostly based on simulations, and focused on wireless ad hoc network scenarios only.

I. INTRODUCTION

The advent of packet radio networks and lead to new concepts such as Mobile Ad hoc Networks (MANETs) or wireless mesh networks. Since the late 1990's, several research efforts have focused on enabling communication in such wireless multi-hop networks in which topology is spontaneous and dynamic. This resulted in different routing protocols being designed for wireless mesh networks and MANET operation, the most prominent to date being OLSR [START_REF] Clausen | RFC 3626, Optimized Link State Routing Protocol (OLSR)[END_REF] and AODV [START_REF] Perkins | RFC 3561, Ad hoc On-Demand Distance Vector (AODV) Routing[END_REF]. OLSR is for instance deployed in many urban wireless community mesh networks around the world [START_REF]Urban Wireless Community Mesh Networks: Freifunk Berlin www.freifunk.net[END_REF].

In practice, the interconnection of wired networks with wireless mesh networks in the Internet gives birth to Autonomous Systems (AS) containing both mesh and fixed networks -such an AS is hereafter denominated compound AS. An approach for routing in a compound AS consists in splitting the AS into several routing domains, and then to use a separate routing protocol for each domain: on one hand MANET-specific protocols in wireless mesh networks, and on the other hand, routing protocols designed for wired networks. An alternative approach consists in using a single routing protocol for the whole compound AS. The former approach requires the presence in the AS of routers with specific hardware and software capabilities, denominated gateways, in order to provide connectivity between different routing domains. The latter approach, explored in this paper, reduces the cost of network maintenance and operation as no such gateways are needed [START_REF] Cordero | Wired/Wireless Compound Networking[END_REF], and naturally avoids bottlenecks created by such gateways. This is, however, at the expense of increasing the complexity of the employed routing protocol, which thus needs to handle the diverse characteristics of wireless mesh and wired networks with the same core mechanisms.

In today's Internet, IGP routing [START_REF] Halabi | Internet Routing Architectures, 2nd Edition[END_REF] is essentially provided by two protocols based on the link state algorithm : OSPF [START_REF] Moy | RFC 2328, OSPF Version 2[END_REF] [10], and ISIS [START_REF] Oran | RFC 1142, OSI IS-IS Intra-domain Routing Protocol[END_REF]. This paper explores the use of OSPF to provide routing in compound ASes, and evaluates its performance on a compound internetwork testbed.

A. Related Work

Leveraging the commonality between OLSR and OSPF (both based on the link state algorithm), operation of OSPF on MANETs has recently been considered. This approach may indeed be cost-effective concerning the integrating of wireless mesh networks in existing infrastructure, as more than 50% of routers in the Internet already run OSPF. Early studies and proposals such as [START_REF] Henderson | A Wireless Interface Type for OSPF[END_REF] [19] paved the way to IETF standardizing several OSPF protocol extensions for operation over ad hoc networks [START_REF] Baccelli | RFC 5449, OSPF Multipoint Relay (MPR) Extension for Ad Hoc Networks[END_REF] [6] [START_REF] Ogier | RFC 5614, Mobile Ad Hoc Network (MANET) Extension of OSPF Using Connected Dominating Set (CDS) Flooding[END_REF]. Various studies have evaluated and compared the performance of these extensions in mobile ad hoc scenarios [START_REF] Baccelli | MPR Techniques with OSPF on Ad Hoc Networks[END_REF], [START_REF] Henderson | Evaluation of OSPF MANET Extensions[END_REF], [START_REF] Baccelli | OSPF over Multi-Hop Ad Hoc Wireless Communications[END_REF], essentially via simulations. Further improvements of these extensions have been proposed in [START_REF] Cordero | MPR+SP -Towards a Unified MPR-based MANET Extension for OSPF[END_REF], [START_REF] Baccelli | Optimization of Critical Data Synchronization via Link Overlay RNG in Mobile Ad Hoc Networks[END_REF], still comparatively evaluated via simulations in MANET only scenarios. In contrast, this paper evaluates OSPF on a real testbed that consists of a wireless mesh network interconnected with a fixed wired network.

II. THE OSPF PROTOCOL AND ITS MANET EXTENSION

OSPF [START_REF] Coltun | RFC 5340, OSPF for IPv6[END_REF], [START_REF] Moy | RFC 2328, OSPF Version 2[END_REF] is a link-state routing protocol for IP networks. This implies that each router maintains a local instance of the Link-State Database (LSDB), representing the full network topology -with the objective of the protocol being that each router should have the same information in its local instance of LSDB and, thus, the exact same view of the network topology. Paths to every possible destination are derived from the Shortest Path Tree (SPT) that every router computes, by way of Dijkstra's algorithm. Although OSPF supports network partitioning in several routing areas, the work described in this paper explores the use of a single area for routing in the whole internetwork.

Routers acquire local topology information and advertise their own presence by periodically exchanging Hello messages with all their 1-hop neighbors (i.e., neighbor sensing).

With such signaling, each router becomes aware of its immediate network topology, i.e. its 2-hop neighborhood. This also allows verification of bidirectional connectivity with 1-hop neighbors (then called bidirectional neighbors).

Each router also explicitly synchronizes its local instance of LSDB with a subset of its bidirectional neighbors. Links between a router and its synchronized neighbors are called adjacencies, and are required to form a network-wide connected backbone, connecting all routers in the network, in order to ensure paths can be computed correctly.

Finally, routers also acquire remote topology information by way of receiving Link State Advertisements (LSA). Each such LSA lists mainly the current adjacencies of the router which generated the LSA. LSAs are disseminated through the entire network in reliable fashion (explicit acknowledgements and retransmissions) via the backbone formed by adjacencies; this operation is called LSA Flooding. Thus, any router which has formed adjacencies must advertise this periodically by way of originating an LSA and performing LSA flooding.

Remote topology information is then used for the construction of the Shortest Path Tree: each router computes the shortest paths over the network graph described in the set of received LSAs it, by way of Dijkstra's algorithm.

A. Packet Types

Routers in OSPF use five types of messages and packets to exchange topology information over the networks, some of them have been already mentioned in this section: Hello packets are used for neighbor sensing, Database Description (DBDesc) packets are exchanged for LSDB synchronization and Link State Advertisements (LSAs) are used for topology reliable flooding and update. After the exchange of DBDesc packets, a router in process of LSDB synchronization may request to its synchronizing neighbor the retransmission of particular LSAs -these requests are sent by way of Link State Request (LSReq) packets. Several LSAs may be sent in a single Link State Update packet (LSU). Several LSA acknowledgements may also be grouped in a single Link State Acknowledgment (LSAck) packet.

B. Interface Types for Wired Links

Rules for flooding and adjacency handling vary for the different interface types supported by OSPF. In broadcast and non-broadcast multiple access (NBMA) interfaces, the flooding procedure is mainly managed by Designated Routers (DRs). A Designated Router is elected among routers whose interfaces are connected to the same link. Such a DR forms adjacencies with all the routers connected to the same link, and it becomes responsible for flooding of LSAs, originated by routers on that link. In point-to-point and pointto-multipoint interfaces, all links are synchronized and all interfaces participate in LSA flooding.

C. MANET Interface Type

The MANET interface type is defined in the extension of OSPF for operation over MANETs. Three different extensions have been standardized by the IETF [START_REF] Roy | RFC 5820, Extensions to OSPF to Support Mobile Ad Hoc Networking[END_REF], [START_REF] Ogier | RFC 5614, Mobile Ad Hoc Network (MANET) Extension of OSPF Using Connected Dominating Set (CDS) Flooding[END_REF], [START_REF] Baccelli | RFC 5449, OSPF Multipoint Relay (MPR) Extension for Ad Hoc Networks[END_REF], each of which specifies mechanisms to optimize topology description, flooding and LSDB synchronization in wireless ad hoc environments.

The experiments carried out used RFC 5449 [START_REF] Baccelli | RFC 5449, OSPF Multipoint Relay (MPR) Extension for Ad Hoc Networks[END_REF]. Wireless interfaces following this specification select a set of Multi-Point Relays (MPRs) among its bidirectional neighbors [START_REF] Clausen | RFC 3626, Optimized Link State Routing Protocol (OLSR)[END_REF]. The set of MPRs selected by the wireless interface of a router must ensure that every packet received from the router reach all 2-hop neighbors of the selecting interface in 2 hops (MPR coverage criterion). A link between a router and one of its MPRs is denominated MPR link.

LSA flooding is then performed through MPRs, meaning that an LSA transmitted (originated or forwarded) by an interface is retransmitted by the MPRs of such interface. Links between interfaces and their MPRs are synchronized and thus become adjacent. Moreover, each interface describes in its LSAs the set of MPRs and MPR selectors. As the set of adjacencies based on MPR selection may not provide a connected subgraph, links from one additional router in the network (denominated synch router) are also declared adjacent to ensure adjacency set connection [START_REF] Cordero | MPR-based Pruning Techniques for Shortest Path Tree Computation[END_REF]. The Shortest Path Tree is then constructed over the set of adjacencies.

III. TESTBED DESCRIPTION

This section describes the characteristics of the employed networking testbed. Section III-A presents the distribution of computers in the testbed and the network topology that they form. Section III-B details the implications of such topology in OSPF routing.

A. Interfaces Configuration and Network Topology

The testbed is composed of 6 fixed computers (routers/hosts) attached to two interconnected networks: a wired network and a wireless network.

1) Physical Topology:

The internetwork connecting these computers was deployed in the Computer Science Lab (Laboratoire d'Informatique, LIX) of École Polytechnique, in Paris (France). Three scenarios -I, II and III-were configured over the resulting internetwork. These scenarios permit to test the communication between computers wless3 and server, for different situations. The physical distribution of computers at LIX is displayed in Figure 1.

Positions of computers do not change, except for the case of wless3, which has a different position for scenario I and for scenarios II and III, as shown in Figure 1. 2) Logical Internetwork Topology: Each scenario corresponds to a specific internetwork topology. Figure 2 indicates the internetwork topology graphs for scenarios I, II and III. In the wired network, computers communicate through the IEEE 802.3 (Ethernet) standard protocol, server is connected with hybrid1 by way of interface eth1 and with hybrid2 by way of interface eth2, as shown in Figure 2. In the wireless network, interfaces communicate through the IEEE 802.11b WLAN standard protocol, and all wireless routers (hybrid1, hybrid2, wless1, wless2 and wless3) use their wireless interface wlan0. The topology that results from wireless reachability among computers hybrid1, hybrid2, wless1, wless2 and wless3 is modified by means of MAC filtering in order to disable links

h 1 ←→ h 2 , w 1,3 ←→ w 2 , w 1,3 ←→ h 2 and w 2 ←→ h 1 .
In scenario III, link h 1 ←→ w 1 is suppressed by disabling interface wlan0 at computer hybrid1.

B. OSPF Routing Configuration

All interfaces use the extended OSPFv3 routing protocol, wired and wireless interfaces using different interface types. Wired interfaces are configured as point-to-point interfaces, as specified in RFCs 2328 [START_REF] Moy | RFC 2328, OSPF Version 2[END_REF] and 5340 [START_REF] Coltun | RFC 5340, OSPF for IPv6[END_REF]. Wireless interfaces are configured as MANET interfaces, as specified in the MPR-OSPF MANET extension (RFC 5449 [START_REF] Baccelli | RFC 5449, OSPF Multipoint Relay (MPR) Extension for Ad Hoc Networks[END_REF]).

1) OSPF Adjacencies and MPRs: According to the specification of OSPF and MPR-OSPF extension, all links in any of the considered topologies for scenarios I, II and III are adjacent. Within the wired network, every point-to-point link is an adjacency. In the wireless network, wireless links are adjacent if they are MPR links. The list of MPRs of every wireless interface, for each scenario, is displayed in Table II.

It can be observed that all links are MPR links, and therefore all are adjacent. In this topology, the presence of a synch router (see section II-C) is thus redundant.

2) OSPF Flooding: Flooding in the wired network is performed through adjacent links, i.e. S ←→ h 1 and S ←→ h 2 . In the wireless network, flooding is performed:

• through the MPR links (from a wireless router towards its MPR), and

• through all links connecting an interface to a hybrid router (hybrid1 and hybrid2).

IV. EXPERIMENTS AND RESULTS

For each scenario (I, II and III), communication between wless3 and server is tested by way of two experiments:

• Transmission of ICMPv6 [START_REF] Conta | RFC 4443, Internet Control Message Protocol (ICMPv6) for the Internet Protocol Version 6 (IPv6) Specification[END_REF] requests (pings) from wless3 to server. The measure of time between the transmission of an ICMP request and its reply corresponds to the Round Trip Time (RTT) of the ping through the evaluated path.

• Transmission of a constant bit rate data UDP flow from wless3 to server. Comparison between packets sent and packets received permits to test the quality of the traversed paths and the wireless links that compose them in each scenario. For a detailed description of these UDP flows, see Table III in the Appendix. Displayed results in this section show the most prominent averaged measurement values over tens of samples (see the Appendix for details); a complete description of the obtained results can be found in [START_REF] Cordero | Compound Wired/Wireless Internetworking with OSPF[END_REF]. The three considered scenarios (I, II and III) are complemented by another scenario in which information is transmitted and measured through the wired link h 1 ←→ S. Results on this scenario are added for completeness. Section IV-A presents the results obtained in both experiments, for each scenario, in terms of quality of wireless links. Section IV-B examines the amount and structure of OSPF control traffic. As expected, the wired link h 1 ←→ S has an almostideal behavior: 100% PDR and no significant delay. The negative impact of wireless links in the path from source to destination is close-to-linear with the number of traversed wireless links, as shown in Figure 3.a: more than 30% of transmitted packets are lost in the first wireless link, and such percentage increases about a 15% per additional wireless link included in the path. Figure 3.b shows that such degradation is also evident in terms of round trip time (RTT). Replies to ICMP requests are immediately delivered through a wired link, but the average and the variation of delays grow with the number of wireless links involved -is in the order of tens of miliseconds for 2 and 3 wireless links.

A. Wireless Mesh Communication

While the impact on communication due to the use of wireless links depends on the specific topology and the network technology that is used, two conclusions can be drawn from these experimental results. As each additional wireless hop in the route of data packets in the network implies a significant degradation of the quality of communication, routing in wireless networks should preserve the principle of shortest (wireless) paths, meaning that the number of wireless links traversed by data packets should be minimized. This confirms the conclusions of simulationbased studies such as [START_REF] Baccelli | OSPF over Multi-Hop Ad Hoc Wireless Communications[END_REF] which highlights the importance of not sacrificing path optimality for less control traffic.

Moreover, in the context of compound internetworks with both wired and wireless links, it is obvious that 'optimal' does not necessarily mean 'the least number of hops', as implicitly used in previous work such as [START_REF] Baccelli | MPR Techniques with OSPF on Ad Hoc Networks[END_REF], [START_REF] Henderson | Evaluation of OSPF MANET Extensions[END_REF], [START_REF] Cordero | MPR+SP -Towards a Unified MPR-based MANET Extension for OSPF[END_REF], [START_REF] Baccelli | Optimization of Critical Data Synchronization via Link Overlay RNG in Mobile Ad Hoc Networks[END_REF]. Indeed, in compound internetworks, it is better for a path to use wired links than wireless links, whenever possible, even if it means more hops in the end. This observation confirms that metrics used should indeed be able to track link quality, and that OSPF sepcifications such as [9] [6] [8] should be completed with a standard way to do that in MANETs. III), and OSPF control traffic. The figures show the structure of such control traffic, in terms of number of bytes, during the first 335 seconds of network operation, i.e., after routers' startup. All interfaces are configured with the same OSPF parameters (see Table IV), in order to facilitate the comparison between control traffic patterns of each of them.

B. OSPF Control Traffic Pattern

1) Hello Packets: The amount of Hello packets sent by each interface is kept constant along the monitored time. As HelloInterval = 2sec, interfaces transmit 2.5 Hello packets per interval of 5 seconds. The length of an average Hello packet is significantly longer in wireless interfaces (Figures 5 and6) than in wired interfaces (Figure 4). For the same number of neighbors, Hellos from server:eth1 have 40 bytes while those from hybrid1:wlan0 have 75.34 bytes (average values). This is due to the fact that Hello packet format in RFC 5449 [START_REF] Baccelli | RFC 5449, OSPF Multipoint Relay (MPR) Extension for Ad Hoc Networks[END_REF] includes additional information about link costs, adjacencies and MPR selection, which is added to the format specified in OSPF [START_REF] Moy | RFC 2328, OSPF Version 2[END_REF] and OSPFv3 [START_REF] Coltun | RFC 5340, OSPF for IPv6[END_REF].

2) LSDB Synchronization: The existence of ongoing LSDB synchronization processes during the monitored time interval can be noticed in the OSPF control traffic structure by way of the presence of Database Description (DBDesc) packets. The fact that these packets are only present, for wired interfaces, in the first part of the monitored interval (from t = 0sec to t = 10sec, as shown in Figure 4 indicates that links become synchronized only when the routers are switched on. In contrast, DBDesc are transmitted in the whole monitored interval for wireless interfaces. This is consistent with the fact that wired links are mostly stable and therefore there is no need to repeat synchronization process after the first LSDB exchange. Wireless links, in contrast, are more prone to packet losses and link failures, and need thus to be synchronized several times during the network lifetime, even in the absence of router mobility. The same phenomenon can be observed with LSRequest packets, which can only be sent during the last phase of the LSDB synchronization process, when the synchronizing neighbors have completed the exchange of DBDesc packets. These observations are also consistent with simulation-based studies such as [START_REF] Baccelli | OSPF over Multi-Hop Ad Hoc Wireless Communications[END_REF] which have studied the impact of link quality on OSPF control traffic. 3) Link State Updates, Requests and Acknowledgements: LSUpdate packets contain one or more Link State Advertise-ments (LSAs). Such LSAs can be either originated by the sending interface, either originated by another interface and flooded (forwarded) by the sending interface. Transmission of LSUpdate packets follows a common pattern in all the interfaces in the internetwork. Thus, pattern consists of periodic peaks followed by time intervals (valleys) in which the number and size of LSUpdate transmissions is lower and roughly constant.

Despite the common pattern in the LSUpdate traffic, several differences can be observed between wired and wireless interfaces. This section concentrates on three particular aspects: peak width, height of valleys between consecutive peaks and transient state (after routers are switched on).

a) Transient state: Immediately after switching on, wireless interfaces transmit a high number of packetsmostly, LSUpdate packets sent in response to LSRequest packets received during the first LSDB synchronization processes in all wireless links (Figures 6 and5, between t = 0sec and t = 50sec). This amount of transmissions involves traffic rates above 220Bps (1.1kB per interval of 5sec), then decreases and stabilizes in a slightly lower level (maximum peak of 130Bps). The opposite behavior is found in wired interfaces (Figure 4), in which the initial transient period of low LSUpdate traffic rate (about 26Bps = 130B 5sec for server:eth1) is followed by a steady period in which the minimum LSUpdate rate is slightly higher (about 29Bps = 145B 5sec for server:eth1). These different behaviors can be explained by the different roles that flooding has over wired and wireless links. Due to their stability, packets sent over wired links are mostly forwarded packets -that is, they come from other interfaces than those involved in the links. In the first instants in which there is no flooding over the network because adjacencies have not been formed in the network and flooding links have not yet been identified, the overall traffic traversing such wired links is temporarily low. The opposite is observed in wireless links.

b) Peak width: Peaks are narrower in wired interfaces (∼ 10sec for server:eth1) than in wireless interfaces (∼ 25sec for hybrid1:wlan0, ∼ 30sec for wless3:wlan0). For interfaces attached to wireless links, there is a high probability that a topology change causes a new topology update before the LSRef resh intervaltherefore, intervals between consecutive transmission of interfaces' topology descriptions are shorter than LSRef resh and the width of the peak increases. In stable wired links, in contrast, intervals between consecutive transmissions are closer to the LSRef resh parameter and, therefore, LSUpdate transmission events are less spread in time.

c) Depth of valleys: Besides the peaks caused by transmission of its own topology description, either periodic or as a reaction to a topology change, two other events may lead an interface to transmit LSAs: (i) forwarding of LSAs originated by other interfaces in the internetwork, and (ii) retransmission of LSAs not acknowledged by their intended destinations. Both additional events explain the presence of valleys with significant traffic rate, i.e., a non-zero minimum level of LSUpdate transmissions in the monitored interfaces. In wired (reliable) links such as server:eth1, such trans-missions are caused by flooding, and involve about 25Bps (127B per interval of 5sec). Wireless interfaces such as wless3:wlan0 have a minimum LSUpdate transmission rate of about 16Bps (80B per interval of 5sec) caused by LSA retransmissions and flooding.

V. CONCLUSION

This paper addresses routing with OSPF across compound Autonomous Systems, i.e., internetworks combining wired networks and wireless mesh networks. Compound internetworks are bound to become a common phenomenon as wireless mesh networks are being deployed and coexist with wired networks. This paper reports on experiments carried out on a compound internetwork testbed confirm prior simulation-based studies analyzing OSPF control traffic over wireless links. We observed that even without mobility and a small number of wireless neighbors, link synchronization involves a continuous and substantial amount of traffic. This observation confirms that reducing the number of synchronized links is a priority when using OSPF over ad hoc networks. Furthermore, data path quality was observed to significantly decrease with each additional wireless hop involved in the path. This observation suggests that wired links should be preferred to wireless links whenever possible, and highlights the need to complement current OSPF specifications with a standard way to track link quality in compound internetworks. The importance of leveraging every possible wired connection indicates that all links in the compound AS should be taking into consideration for route computation. As this implies a large number of gateways between routing domains when several routing protocols are used in the AS, and this is costly, a single routing protocol is prefered in order to avoid gateways.

APPENDIX

• Testbed Hardware -Wired interfaces: DECchip 21140.

-Wireless interfaces: Broadcom BCM4306 WLAN. -RTT results averaged over 60 iterations.

-ICMPv6 requests did not overlap with UDP flows.

• Setup for Control Traffic Measurements -Routers were switched within [0, 2]sec.

-Results averaged over 84 iterations.

-ospf6d daemon restarted in each iteration.

-UDP flows: see Table III.

• OSPF Parameters: see Table IV

Fig. 1 .

 1 Fig. 1. Computers position over the plan of LIX.

Fig. 2 .

 2 Fig. 2. Considered topologies for scenarios I, II and III.

Figures 3 .

 3 Figures 3.a and 3.b display the results of the performed experiments, in particular the delay for ICMP requests (pings) and the packet delivery ratio of CBR UDP data flows. Figures 3.a and 3.b indicate the degradation of the quality of communication between routers wless3 and server as the number of wireless links between them increases.As expected, the wired link h 1 ←→ S has an almostideal behavior: 100% PDR and no significant delay. The negative impact of wireless links in the path from source to destination is close-to-linear with the number of traversed

Fig. 3 .

 3 Fig. 3. PDR of UDP flows and RTT of ICMP requests, depending on the number of wless hops.

Fig. 4 .

 4 Fig. 4. Control traffic overhead at server:eth1.

Figures 4 ,

 4 Figures 4, 5 and 6 display the evolution of OSPF control traffic transmitted by wireless interfaces wless3:wlan0 and hybrid1:wlan0, on one side, and the wired interface server:eth1, on the other. The five packet formats used in OSPF (Hello, LSUpdate, LSRequest, LSAck and DBDesc, see section II-A) can be distinguished in these figures. Measurements were taken with the topology of scenario I, each point indicating the number of bytes sent within an interval of 5 seconds. The traffic load of the internetwork was composed of a CBR UDP data traffic flow from wless3 towards server (see TableIII), and OSPF control traffic. The figures show the structure of such control traffic, in terms of number of bytes, during the first 335 seconds of network operation, i.e., after routers' startup. All interfaces are configured with the same OSPF parameters (see TableIV), in order to facilitate the comparison between control traffic patterns of each of them.1) Hello Packets: The amount of Hello packets sent by each interface is kept constant along the monitored time. As HelloInterval = 2sec, interfaces transmit 2.5 Hello packets per interval of 5 seconds. The length of an average Hello packet is significantly longer in wireless interfaces (Figures5 and 6) than in wired interfaces (Figure4). For the same number of neighbors, Hellos from server:eth1 have 40 bytes while those from hybrid1:wlan0 have 75.34 bytes (average values). This is due to the fact that Hello packet format in RFC 5449[START_REF] Baccelli | RFC 5449, OSPF Multipoint Relay (MPR) Extension for Ad Hoc Networks[END_REF] includes additional information about link costs, adjacencies and MPR selection, which is added to the format specified in OSPF[START_REF] Moy | RFC 2328, OSPF Version 2[END_REF] and OSPFv3[START_REF] Coltun | RFC 5340, OSPF for IPv6[END_REF].2) LSDB Synchronization: The existence of ongoing LSDB synchronization processes during the monitored time interval can be noticed in the OSPF control traffic structure

Fig. 5 .

 5 Fig. 5. Control traffic overhead at wless3:wlan0.

Fig. 6 .

 6 Fig. 6. Control traffic overhead at hybrid1:wlan0.

  Table I indicates the network interfaces of each computer. For more details about computers' hardware, see the Appendix.

	Computer	Abbr. Wired ifs. Wireless ifs.
	server	S	eth1, eth2	-
	hybrid1	h1	eth1	wlan0
	hybrid2	h2	eth1	wlan0
	wless1	w1	-	wlan0
	wless2	w2	-	wlan0
	wless3	w3	-	wlan0

TABLE I NETWORK

 I INTERFACES OF TESTBED COMPUTERS.

TABLE II MPRS

 II SELECTED BY EACH WIRELESS INTERFACE, FOR EACH SCENARIO.

  Operating System: Ubuntu v.10.04 with kernel 2.6.32. -Routing Protocol Implementation: ospf6d daemon of Quagga/Zebra routing suite v.0.99.15. * Wired interfaces: Point-to-point. * Wireless interfaces: MANET (RFC 5449 [9]). • Setup for PDR and RTT Measurements -Routers were switched on within [0, 2]sec. -PDR results averaged over 84 iterations. -ospf6d daemon NOT restarted in each iteration.

	-UDP flows, started 60sec after ospf6d daemon switch-
	on:
	Nominal sender bit rate	100 pkts/s
	Packet payload	1024 bytes
	CBR real traffic rate	∼ 300 kbps
	Flow duration	5 min/flow

• Testbed Software -

TABLE III CHARACTERISTICS

 III OF TRANSMITTED UDP FLOWS.

  .

	Name	Value
	HelloInterval	2 sec
	DeadInterval	10 sec
	RxmtInterval	5 sec
	AckInterval	2 sec
	MinLSInterval	5 sec
	MinLSArrival	1 sec
	LSRefreshInterval 60 sec
	TABLE IV
	OSPF PARAMETERS.


M Tournier

L Reveret

Principal Geodesic Dynamics

Keywords: (according to ACM CCS): I, 3, 7 [Computer Graphics]: 3D Graphics and Realism-Animation

Figure 1: One foot balance. 3 control features are present here: right foot, Center of Mass (CoM), and left hand (red). The system automatically adjusts the left leg position to maintain the CoM above the right foot while the user manipulates the left hand position.

Introduction

While human body exhibits numerous degrees of freedom for producing motion, bio-mechanical studies showed they tend to be actuated in a highly-correlated way, resulting in similarly coordinated body motion. This fact has been typically exploited by recent successful approaches based on dimension reduction of human motion. This redundancy offers an opportunity to reduce the total complexity of animating a virtual character, by capturing this structure in a simple parametric model. It allows to generate better-looking animations at a lower computational cost. We propose in this paper an approach to use these reduced parameters directly at the core of a physically-based animation context by using them as generalized coordinates.

In the context of character animation, recent works have shown the potential of mixing dimension reduction with physical simulation. What we propose here is a complete derivation of their integration into a coherent framework working in real-time and with stable control. Our physical simulation includes constraint-based unilateral ground contact model. Constraint-based model is favored over the penalty-based approach to avoid the risk of instability due too strong stiffness. Our approach allows to implement agile change of foot support during simulation, on user request, without planning the contact location explicitly. In order to maintain a real-time simulation while preserving important dynamic invariants, such as angular momentum, we used an explicit integrator based on variational geometric integrators.

To summarize, our contributions are the followings:

• an explicit time integrator for a reduced dimension pose model based on variational geometric integrators for articulated bodies, • a simple quadratic programming control framework, providing trade-off between physical-simulation and kinematics control using a metric interpolation and real-time animation.

The remainder of the paper is organized as follow:

Section 2 reviews related works in the domain of dimension reduction in the context of both machine learning approach and physically-based animation. Recent works in the domain of task-space control are also discussed. Section 3 presents the algorithm for learning and approximating the pose manifold of a motion capture. Section 4 applies these ideas in the context of physically-based animation of characters using variational integrators. Section 5 on control develops the physical modeling to implement motion control. Examples of controllers for features tracking are presented in 6, leading to results for a balance controller. Finally, benefits and limitations are discussed in section 7.

Related Work

Dimension Reduction in Machine Learning

Recent advances in the machine learning theory allowed to design systems that capture geometric features as well as modeling time evolution on existing motion capture data. They allow synthesizing plausible motion while still retaining good control.

The idea of using machine learning to compute a higherlevel parametrization of a motion given real examples has been widely explored over the past decade. [RCB98] used Radial Basis Functions (RBF) to learn an interpolation space from examples. [START_REF] Mukai | Geostatistical motion interpolation[END_REF] proposed to use geo-statistics to perform better interpolation in a control space. Brand and Hertzmann [BH00] propose a cross-entropy optimization framework to learn HMM for both structure and style from motion capture data. [START_REF] Grochow | Style-based inverse kinematics[END_REF] cast the Inverse Kinematics (IK) problem as the optimization of the likelihood of a Probability Density Function (PDF) over poses, knowing constraints. They propose to learn this PDF using a Gaussian Process Latent Variable Model, a generalization of PCA and RBF models that is computed on motion features (e.g. angles, velocities). Then, the corresponding pose likelihood function is optimized under geometric constraints.

All these methods exploit ever more complex non-linear statistical techniques to represent a statistical distribution over poses, possibly using regression against a feature space or directly optimizing the likelihood in which case the computational cost is high.

Dimension Reduction in Physically-based Animation

The benefits of dimension reduction have also been exploited in the context of physically-based animation, since they allow to reduce computational requirements while preserving most dynamic features of the full-dimension simulation.

[PW89] pioneered the use of modal analysis in the context of computer graphics. The basic idea of this approach is to decompose a mechanical system with numerous, coupled mechanical degrees of freedom into an equivalent set of mechanically-independent one-dimensional degrees of freedom, usually called modes. The set of all modes forms the modal basis which can be seen as an alternate representation of the system DOFs.

Due to its numerous advantages, modal analysis has spawned a wide variety of research works ( [JBP06, TLP06, WST09], among others) building on this basic principle. [START_REF] Kry | Modal locomotion: animating virtual characters with natural vibrations[END_REF] proposed an extension of this technique to articulated rigid bodies: rather than performing the modal analysis on a mesh-supported, linearized FEM, the authors applied the modal decomposition to the dynamics of an articulated rigid body around a rest pose, obtaining new mechanicallyindependent angular degrees of freedom around a rest pose. This was the basis for animating locomotion behaviors as combinations of natural vibration modes. Simulation including contacts have been recently proposed in [JL11] and [NCNV * 12]. However, in this case, modal analysis has the drawback that the joint stiffness must be known for the modal analysis to make sense. Instead, the reduced basis we obtained is not dependent on hand-chosen stiffness/damping parameters, but only on motion capture data. Besides, expressing the dynamics in our reduced pose parametrization is not restricted to small displacements.

Closer to data-driven statistical approach, [SHP04], and more recently [START_REF] Chai | Physically valid statistical models for human motion generation[END_REF], present frameworks where statistical approaches and physically-based animation are integrated. Motion is parameterized with usual kinematics degrees of freedom. In their case, motion manifold learned from PCA acts as a constraint to plausible poses. Such an approach implies an optimization over a period time of several frames. We are proposing a different approach where latent variables induced by manifold learning are directly used as generalized coordinates avoiding a costly optimization.

[YL08] present an alternate physically-based approach featuring dimension reduction for character animation. The goal of this work is to construct a data-driven basis of torques from motion capture data, then extract the coordinates corresponding to the least actuated torques. This near-unactuated basis is later used to synthesize upper-body, physicallybased perturbations on top of existing motion capture animations. The main assumption is that the unactuated coordinates tend to be much more compliant, in the case of external disturbances, than the actuated ones. While this technique shows good results, its drawback compared to expressing the dynamics in the reduced dimension pose space is that it requires inverse dynamics, which is challenging to setup in the case of contacts, probably explaining that the motion perturbation is limited to the upper-body in their context.

Motion Control

The goal of this paper is not to propose new control algorithms, such as ones implementing walking. Our goal is rather to show how control can be implemented using our new formulation of the dynamics of articulated characters. We follow a task-space control approach as described in [AdSP07,MZS09,JYL09,LMH10,MLH10]. In this formulation, kinematic goals are formulated in some abstract task, or feature space, and the system solver is responsible for optimizing these goals using physics laws as a constraint at each time-step. This approach can be seen as a synthesis of the advantages of joint-space control and space-time optimization while mitigating their drawbacks: instead of solving a large, global optimization problem preventing interactivity, a smaller, easier problem is solved at each time-step. Besides, the burden of choosing and adapting individual control gains is shifted from the joint space to an abstract task-space.

In [LMH10, MLH10], the kinematic goals are expressed on accelerations, by forming a quadratic objective function. This objective is optimized under the dynamic equations and approximate contact constraints (no complementarity). A feature priority solving strategy is proposed, preventing certain control objectives to compete with each others. While this approach is very general and well-grounded for character control, it could arguably be improved especially by the use of a velocity/impulse formulation of dynamics and features, removing the need for second derivative which can be difficult to compute. Last, but not least, a reduced dimension motion model could improve the overall efficiency of this method significantly. These are the points we will address in the remaining of this paper.

Reduced Manifold for Kinematics

The configuration space of an articulated character, parameterized by its joint orientations, is not an Euclidean vector space. This is the key limitation which prevents from applying safely standard analysis such as PCA. However, the pose manifold SO(3) n possesses a natural Riemannian structure, for which generalizations of Euclidean algorithms have been proposed, including multi-resolution analysis, interpolation, and statistics. Building on this property, [FLPJ04] proposed an extension of PCA to the case of Riemannian manifolds, first in the case of Lie groups having a compatible Riemannian structure, then on symmetric Riemannian manifolds.

Following [TWC * 09] in the context of character animation, we have found the PGA algorithm to be a natural, simple dimension reduction technique, suitable for analyzing rotational pose data. We quickly recall main results in the context of character animation in order to introduce notations.

Geodesics on Lie Groups

Geodesics on a smooth manifold arise from a Riemannian metric, a smooth operator computing the norm of tangent vectors and the length of curves: geodesics are locally length-minimizing curves. Lie groups, on the other hand, are smooth manifolds where the group operations are smooth. A good introduction to these topics in the context of robot kinematics can be found in [MLS94].

In the case of a Lie group with a compatible Riemannian metric such as SO(3), the metric and algebraic exponential coincide. Given such a Lie group G, the geodesic starting at a ∈ G with initial tangent vector V ∈ g (in body coordinates) is given by γ a,V = a. exp(αV ), for α ∈ R (cf. Figure 2). The geodesic distance between two points a and b of G is then given by:

dist(a, b) = log(b -1 a)
In terms of Lie group operations, the projection of x ∈ G on a geodesic starting at a ∈ G with body tangent vector V ∈ g is thus:

π a,V (x) = a exp(α V ) (1) with α = argmin α∈R dist (x, a exp (αV ))

Principal Geodesic Analysis

We now quickly outline the PGA algorithm for orientation data. More in-depth discussions can be found in [FLPJ04, SLHN10] for the general case. In essence, the PGA algorithm generalizes PCA using the geodesic distance and projection instead of the Euclidean ones. Given orientation samples, the first step is to compute their intrinsic mean µ ∈ SO(3), which minimizes the sum of geodesic distances to the samples [START_REF] Moakher | Means and averaging in the group of rotations[END_REF]Pen06]. The intrinsic mean provides a coordinate-invariant description of the data average. In contrast, an Euclidean mean of Euler angles or exponential maps depend on a particular choice of reference coordinate frame. Next, a set of mutually orthogonal geodesic directions is computed by maximizing the geodesic projections of the samples onto the corresponding geodesics.

Dimension reduction is obtained by selecting the k ∈ N most significant geodesic directions, where k is usually determined using a variance criterion. As the PGA algorithm only makes use of Lie group operations, its extension to the case of pose data, i.e. the product group SO(3) n with n ≥ 1, is straightforward.

An approximate computation of the geodesic directions can be obtained by applying a PCA over the log-linearized samples at the intrinsic mean µ. While this approach somehow mitigates the advantages of PGA, we found it more practical to use and sufficient for animation purposes, while retaining the coordinate-invariance of the analysis. The differences between linearized and exact PGA are discussed in detail in [SLHN10].

We now describe how to apply PGA on motion capture data in order to construct reduced dimension skeleton kinematics.

Reduced Kinematics

We begin with an existing motion capture sequence represented as a set of m skeleton configurations (g i ) i≤m ∈ G sampled over time. In practice, we found that breakdance sequences provide interesting variability in poses, with an associated set of 10 to 20 geodesics offering a range of motion large enough to allow generalization to other motion. Each configuration is a pair (c, p) ∈ G = SE(3) × SO(3) n , where n ∈ N is the number of joints, c is the rigid transformation describing the configuration of the root bone, and p contains the relative orientation of each joint with respect to its parent.

We will denote the pose space of relative joint orientations by P = SO(3) n . We apply PGA on the pose samples (p i ) i≤m to obtain the intrinsic mean µ ∈ P, and a set of k ∈ N geodesic directions V j ∈ p describing the pose submanifold. With this decomposition, the new pose DOFs are k scalars (α j ) j≤k ∈ R describing coordinates over the k principal geodesics. Practical computation of the intrinsic mean µ and geodesic directions V j are detailed in [TWC * 09].

We define the reduced pose parametrization as a smooth function r : R k → P that reconstructs a pose given k geodesic coordinates, α = (α j ) j≤k . Two mappings may be chosen for pose reconstruction, known as the canonical coordinates of the first and second kind [MLS94]. In practice, both have provided similar results for our purpose. We thus restrict to the first one:

r(α) = µ exp k ∑ j=1 α j .V j
We obtain the complete character forward kinematics, including the root DOFs, by composing r with the usual forward kinematics for articulated rigid bodies. We summarize this mapping as the following function:

f : SO(3) × R 3 × R k → SE(3) n+1 (3)
where the configuration space G := SO(3)×R 3 ×R k is a Lie group describing respectively the absolute orientation, absolute position and reduced pose coordinates of the character. f computes the absolute configuration of each of the k + 1 bones of the skeleton.

The differential of f is obtained by the chain rule. Note, however, that the derivative of the exponential does not assume a simple formula as it does for the real scalar case, due to the non-commutativity of the rotation group. See [MLS94] for appropriate derivation.

Principal Geodesic Dynamics

Having described how to use a PGA-based reduced pose model in a kinematic animation context, we now move on to the physically-based animation of a character, using this reduced model as generalized coordinates. In this section we motivate and describe a time-stepping scheme for physically animating a virtual character parameterized using PGA. We derive a geometric integrator based on [START_REF] Kobilarov | Lie group integrators for animation and control of vehicles[END_REF]. Finally, we propose an approximated, explicit time-integrator, for use in a real-time simulation.

The use of reduced coordinates is natural in our context since we already have the reduced coordinates pose mapping as the PGA pose parametrization. In this case the constraint forces are implicit, and no constraint drift can occur by design. However, this formalism imposes the derivation of dedicated equations of motion, and notably to compute the reduced mass tensor and Coriolis forces according to the reduced coordinates mapping. In this context, the lowdimension of the PGA parametrization will be an advantage since it will keep the computational cost reasonable while still enforcing natural poses.

Geometric Integrator Derivation

We now derive the discrete equations of motion for a PGAparametrized virtual character. For this, we use the geometric, variational integrator methodology presented in [START_REF] Kobilarov | Lie group integrators for animation and control of vehicles[END_REF]. This work addressed the case of a single rigid body and consequently made the assumption of the leftinvariance † of the Lagrangian. Instead, we derive equations for the general case of multiple articulated rigid bodies, where such an assumption does not hold. Unless stated otherwise, all the tangent maps and vectors used in the following are expressed in body coordinates, as opposed to spatial coordinates [MLS94].

Let f : G → SE(3) n+1 be the forward kinematics as defined in Equation (3), and J(g) its Jacobian matrix. We assume that the inertia tensors for the n + 1 bones composing the skeleton are given as a positive definite, block-diagonal matrix M of dimension 6(n+1). The kinetic energy T is then defined as:

T (g, v) = 1 2 v T J(g) T M J(g)v
where g ∈ G is the configuration, and the velocity v ∈ g belongs to the Lie algebra. The above expression defines the generalized mass matrix, or inertia tensor as:

M(g) = J(g) T M J(g)
The potential energy V (g) is simply the gravitation potential.

Variational Integrator Assuming a Lagrangian L = T -V is given, the Hamilton equations of motion can be derived from the Hamilton-Pontryagin (HP) principle [KYT * 06]. By properly discretizing the HP principle through a quadrature, one can obtain a symplectic integrator of arbitrary accuracy (called variational integrator) with excellent momentum and energy conservation properties. In the case of Lie groups, [START_REF] Kobilarov | Lie group integrators for animation and control of vehicles[END_REF] propose a first-order accurate quadrature of the HP principle in terms of a group difference map τ : G → g:

δ T 0 L(g, v) + p T ( ġ -v) dt ≈ (4) δ N ∑ k=0 hL (g k , v k+1 ) + p T k+1 τ g -1 k g k+1 -h v k+1 = 0
where the momentum p T ∈ g * belongs to the Lie coalgebra, and h is the time step. In practice, we used the logarithm as the group difference map. The stationary conditions described by the above discretized HP principle result in a symplectic update rule for g, v and p, summarized in the following non-linear system (see Annexe A for more details):

∂L ∂v (g k , v k+1 ) dτ k+1 = p T k dτ k Ad d k + h ∂L ∂g (g k , v k+1 ) (5a) g k+1 = g k τ -1 (h v k+1 ) (5b) p T k+1 = ∂L ∂v (g k , v k+1 ) (5c)
where d k+1 = g -1 k g k+1 , dτ k := dτ(d -1 k ) and Ad is the † i.e. constant in body coordinates adjoint ‡ representation of G. Forcing can be incorporated through an extended variational principle [START_REF] Kobilarov | Lie group integrators for animation and control of vehicles[END_REF]. If f k is an external force applied at the k th time step, the integrated force h f k is added to the velocity update (5a):

∂L ∂v (g k , v k+1 ) dτ k+1 = p T k dτ k Ad d k + (6) h ∂L ∂g (g k , v k+1 ) + f T k
Equations (5b) and (5c) are left unchanged.

Approximations for Real-Time Simulation

Two non-linear expressions in v k+1 complicate the integration update in equation ( 6). We propose two approximations to turn the update into a linear system, at the expense of symplecticity. Our experimental results show however that angular momentum is nearly conserved, which is not the case when using a non-Lie group integrator.

Quadratic Forces

The presence of the quadratic forces [MLS94] ∂L ∂g (g k , v k+1 ), containing Coriolis and centrifugal forces, is problematic since this term is quadratic in v k+1 . To avoid this issue, we approximate it by using the velocities from the previous time step, thus making it completely explicit:

∂L ∂g (g k , v k+1 ) ≈ ∂L ∂g (g k , v k )
In our simulator, we compute this term using central finite differences.

Difference Map Derivative

The second cause of non-linearity in the above variational update is the presence of dτ k+1 in the left-hand side of Equation (6), and this term depends non-linearly on v k+1 . As above, we approximate this term by reusing the velocity from the previous time step, which corresponds to d k+1 ≈ d k . Under these approximations, the velocity update becomes the following linear system:

∂L ∂v (g k , v k+1 ) = p T k dτ k Ad d k + (7) h ∂L ∂g (g k , v k ) + f T k dτ -1 k
Note that if a more accurate predictor is available for v k+1 , the two approximations presented above can be adapted in a straightforward manner. ‡ For a Lie group G, the adjoint Adx : g → g maps body coordinates to spatial coordinates corresponding to a tangent vector at x ∈ G.

Evaluation

We now evaluate the proposed integrator. To show the advantages of Lie group integration, we compare our integrator with one of the following form:

∂L ∂v (g k , v k+1 ) = p T k + h ∂L ∂g (g k , v k ) + f T k (8)
Equation ( 8) is a naive adaptation of the Euclidean symplectic Euler integrator to the Lie group case, with explicit quadratic forces for the sake of comparison. The extra terms dτ k Ad d k and dτ -1 k in Equation ( 7) can be seen as appropriate change of frames for taking the Lie group curvature into account when summing external forces and previous momentum. A graphical interpretation of this can be found in [START_REF] Kobilarov | Lie group integrators for animation and control of vehicles[END_REF].

We compare the behavior of these two integrators on the following scenario: the character is in a gravity-less environment, we apply a short impulse on its right hand, along the positive Z axis. We choose a quite large time-step for an explicit integrator: h = 0.1s.

The accompanying video shows the trajectory of the Angular Momentum (AM) vector at the center of mass position. The Euclidean integrator clearly shows instability that our integrator does not exhibit, as seen in Figure 3. Another experiment is performed with a fall under gravity onto an infinite inclined slope. Again, the conservation of momentum appeared to be preserved by our integrator while some artifacts arose for the Euclidean integrator.

Figure 3: The Angular Momentum (blue) is displayed after a short lateral impulse in a gravity-less environment. Our Lie group integrator (left) features a much better AM conservation than a non Lie-group one (right).

Control Formulation

Now that a physically-based animation model has been derived, we turn to the problem of motion control.

Quadratic Programming

We describe a character control framework in the same spirit of [LMH10], but at the velocity/impulse level instead of force/acceleration. Equation ( 7) performs an explicit velocity update as a linear system, which can be summarized in a compact velocity-impulse formulation Mv = p. If we incorporate contact impulses, non-penetration and complementarity constraints into this system, we obtain the Karush-Kuhn-Tucker (KKT) conditions of a convex Quadratic Program (QP). We now formulate motion control by altering this QP.

The constrained dynamics of a character can be summarized as the following QP:

v k+1 = argmin Av≥b 1 2 v T Mv -p T v ( 9
)
where v is the velocity coordinates, p is the net momentum (i.e. the left-hand side in Equation ( 7)), M is the mass tensor, and A, b describe unilateral constraint geometry. We modify the original quadratic form, using a control quadratic form 1 2 v T Qv + c T v, in order to obtain a different behavior:

v k+1 = argmin Av≥b 1 2 v T Mv -p T v dynamics + 1 2 v T Qv + c T v control ( 10
)
The control form may represent any soft kinematic constraint one wishes to satisfy on velocities, expressed as a quadratic form. We will see how this form is constructed in the next paragraphs. To enrich the user control, a simple interpolation parameter can be employed for the two quadratic terms to balance between a purely physical simulation, or better follow the kinematic constraints at the expense of physical realism.

One should note that modifying the original quadratic form in the above way effectively amounts to adding implicit forces of the form Qv + c, as can be seen on the corresponding KKT conditions. If the control form involves the root DOFs of the character, the associated forces directly act on the character root.

We now describe how to construct the control form in terms of control features.

Features Tracking Formulation

In practice, the control form is the aggregation of quadratic error terms, each one controlling a specific velocitydependent feature. Let F R d be an abstract feature space of dimension d. We define a feature γ as an affine map γ:

γ : g → F γ(v) = Γv -γ
We call Γ ∈ M d,n the feature matrix, and γ ∈ F the feature desired velocity. Given a set of m ∈ N features (γ i ) i≤m , we build a quadratic form by simply summing all the feature squared norms:

1 2 v T Qv + c T v := 1 2 m ∑ i=1 ||γ i (v)|| 2 More precisely, Q = ∑ i Γ T i Γ i , c = -∑ i Γ T i γ i
, and unneeded constant terms are dropped. Features can be given more importance in the final control form by weighting them accordingly. Such a weight is implicit in the definition of each feature, in the above formula.

We now illustrate the above formulation on a simple IK example. Let a function e : G → R 3 describe the error between an end-effector position and its desired position. The error for the current configuration g k is e k := e(g k ). We obtain an estimation of e k+1 based on a first-order approximation at the current time step:

e k+1 ≈ e k + h Je k v
where Je k is the Jacobian matrix of e at g k . We would like e k+1 to be as small as possible, thus we try to enforce:

Je k v = - e k h
The corresponding control feature is described by the pair (Je k , e k /h). In this example, we ask the error to be zero at the next time step, which corresponds to a proportional gain of φp = 1/h. Smaller gains can also be provided, as well as derivative gains, in order to obtain a general Proportional-Derivative (PD) model for control features:

Je k v = -φp e k -φ d ėk ( 11
)
where ėk = Je k v k is the error derivative for the current time step. As one would expect, the derivative gain φ d tends to damp oscillations around the control target, at the expense of convergence speed. In the case where the error function is defined on a Lie group (e.g. when controlling orientation), the Lie group PD control described in [BM95] can be used instead. An example of the IK control feature is shown in Figure 4.

We now present some control results obtained using our system.

Examples of Controllers

Figure 4: Simple IK control feature. The desired position (red) is interactively manipulated by the user while the character maintains balance.

Looking at Targets

Due to the full-body correlations captured by the PGA pose model, we have found the head rotational motion to be sometimes too important and thus penalizing for visual quality.

We thus implemented the following target tracking feature in order to gain an intuitive control on the head orientation.

The absolute head orientation is defined by a mapping q : G → SO(3). We compute the desired orientation as a direct orthogonal basis of R 3 as follows:

• The first basis vector b 1 is the normalized view vector w ∈ R 3 , constructed as the difference between the view target and the head center • The second basis vector b 2 is the projection of the world up vector u ∈ R 3 onto the plane orthogonal to b 1 , normalized • The third basis vector b 3 is obtained by taking the crossproduct of the two first vectors, finalizing the basis Such construction is obviously ill-defined when the view and up vectors coincide. Should this situation happen, we simply skip this feature in the final control form.

Once the desired orientation q ∈ SO(3) is computed, we define an error function e(g) = log(q -1 q(g)) ∈ R 3 and obtain the corresponding control feature as described in 5.2. Figure 5 shows an example use of this feature control.

Balance Control

Let us now turn to balance-related feature control. As shown in [MZS09], the linear and angular momenta are two highlevel quantities of particular interest in this context.

Center of Mass (CoM)

The CoM position is given by a function c : G → R 3 defined as:

c(g) = n+1 ∑ i=1 m i c i (g)
where n + 1 ∈ N is the number of rigid bodies in the skeleton, c i computes the position of the CoM for the i th body, and m i is the mass for the i th body.

Given a target position for c, an error function is easily derived from the above definition. Again, we obtain the corresponding velocity feature as described in 5.2.

Angular Momentum (AM)

In a feature-based control framework, it is possible to directly control the AM [LMH10], as opposed to an indirect AM control through the Zero Momentum Point [MZS09]. The AM with respect to the CoM c is defined as:

a = n+1 ∑ i=1 R i I i ω i + n+1 ∑ i=1 m i (c i -c) × ċi ∈ R 3
where R i ∈ SO(3) is the absolute orientation of the i th body, I i is the body-fixed inertia tensor of the i th bone, ω i is the i th body-fixed angular velocity. The AM is linear in the rigid body velocities, which are themselves linear in the generalized velocities. Thus, the AM may be rewritten as a linear function of the generalized velocities:

a = H v
Thus, H is used as the AM feature matrix. We set the desired AM value to zero in order to prevent the character from undergoing destabilizing rotational motion.

Support Center

Instead of fully controlling the CoM position, it is usually sufficient to keep its projection on the ground inside the support polygon [JYL09, LMH10]. To do so, we record the position of the character contacts with the ground at each time step. We define the center of support s ∈ R 3 as the mean of these contact points (cf. Figure 8). A control feature is added so that the CoM projection on the ground matches the center of support.

Results

Our balance controller uses the following feature set:

• Head and CoM projections should be at the support center, • Feet should stay on the ground, • Angular momentum should be minimal, • CoM height.

We found the head projection criteria to be a simple yet effective way to maintain the character upright. Figure 1

Discussion and Limitations

Efficiency As expected, the dimension reduction offered by the PGA allows to improve efficiency, since all the examples we proposed run in full real-time, at around 100Hz on a quad-core machine. We believe these performances to be improvable by engineering a better optimized code. Joint Limits Though not described here for clarity, we have also added handling of joint limits. Our approach is based on a sampling of values observed in the learning joint orientation data. A Minimum Volume Enclosing Ellipsoid is fitted to these data as a compact, analytical description of the allowed joint space. Unilateral constraints preventing the joint orientation from leaving this ellipsoid were incorporated in the control framework described in section 5. A conceptually similar, but more complex approach can be found in [HUH02].

Robustness Though the balance controllers we present performed generally well under user interaction, summing control terms in an unprioritized way will necessarily be subject to robustness problem in the case of a user insisting too much. We did not implement yet a prioritized optimization strategy as found in [LMH10], though it could have largely benefited from our reduced dimension strategy. This kind of framework produces much more robust controllers than simply summing objective terms, as the core balance objectives can not be perturbed by user IK requests. This could also improve the robustness of our controllers to impulsive perturbations.

Generalization Principal geodesic DOFs are completely dependent on the choice of the motion capture sequence used for learning. It thus questions the generalization capacity of such an approach. In our experiments, we mainly used the breakdance sequences for learning as shown in the accompanying video. It shows good properties for the interaction scenario we provide. The computation of PGA is extremely fast (less than a second in the example shown). The learning phase is not a limitation and different learning data could be tested on the fly if the resulting geodesics seem not appropriated to the type of motion edition wished by the user. More complex behavior such as locomotion is beyond the intended scope of this paper. As shown by numerous previous studies, simulation of locomotion requires careful controllers (typically finite state machines are needed). Our goal here was to firstly show that features-based controllers could be implemented in our Principal Geodesic Dynamics framework. The development of more complex controllers implementing locomotion is a natural extension of our work.

Conclusion

A Quadratic Programming control approach, exploiting a PGA reduced pose parametrization, has been presented. This approach allows to easily add small dynamic effects to the kinematic manipulation of a virtual character. We have experienced the system to be stable enough even when the contacts scenario is not known in advance. While such a control framework presented above can produce good visual results, it is not totally physically correct since the character is externally actuated through the features tracking. A more satisfying, but also more difficult approach, is to only use the internal actuators of the character to achieve a given task expressed using velocity features. Doing this is much more difficult from a theoretical point of view, as the true task-space control problem under unilateral constraints is non-convex, and thus much harder to solve. Some algorithms have been proposed for solving such problems, such as a Sequential Quadratic Program. We did not try these methods yet, looking for real-time solutions in the first place. However, although these iterative algorithms are computationally-expensive, our reduced pose model may finally provide a low enough dimension reduction to make their real-time use possible.

Annexe A: Discrete Geometric Integrator

Let d k+1 = g -1 k g k+1 . We now derive the stationary conditions for momentum, velocity and position in Equation (4).

Momentum Taking variations in p gives:

N ∑ k=0 δp T k+1 (τ(d k+1 ) -h v k+1 ) = 0
The fundamental lemma of the calculus of variations implies that for all k ≤ N:

δp : τ(d k+1 ) = h v k+1 ∈ g (12)
Velocity Taking variations in v gives:

N ∑ k=0 ∂L ∂v (g k , v k+1 ) -p T k+1 δv k+1 = 0
The same argument implies that for all k ≤ N:

δv : ∂L ∂v (g k , v k+1 ) = p T k+1 ∈ g * (13)
Position Let us first remark that τ(x -1 ) = -τ(x), and let dτ k := dτ(d -1 k ). Taking variations in g with fixed end-points (δg 0 = δg N+1 = 0) gives:

Figure 2 :

 2 Figure 2: The geodesic projection of point x over the geodesic exp a (αV ) (blue) corresponds to the point where the geodesic segment in green has minimal length.

Figure 5 :

 5 Figure 5: An example of head orientation control: the character is looking at the red target.

Figures 6

 6 Figures 6 and 7 compare the resulting behaviors when the AM control is turned off and on.

Figure 6 :

 6 Figure 6: One foot balance, without AM control: the character makes broad moves and quickly falls down. The red marker shows the desired CoM position.

Figure 7 :

 7 Figure 7: One foot balance, with AM control: the character adjusts his arms and moves in a balanced way. The red marker shows the desired CoM position.

Figure 8 :

 8 Figure 8: The contact points (red) are used to define the support center (blue) as their mean. The CoM projection (green) is controlled by setting its desired value to the blue point.

  and Figure 9 show examples of one-foot balance control. The accompanying video shows extended examples of interaction, with on-line change of support foot controlled by the user. This behavior is simply obtained by biasing the CoM and head projection targets toward the supporting foot position. The overall posture changes (standing, crouching, jumps) of the character are implemented through change in target height of the CoM.

Figure 9 :

 9 Figure 9: One foot balance, while tracking an IK target.

  , v k+1 ) δg kp T k+1 dτ k+1 δd -1 k+1 = 0 (14)Noting that δd -1 k+1 = -Ad d k+1 δg k+1 + δg k , we split the sum by grouping terms in δg k and δg k+1 together. After renumbering and exploiting the fixed endpoints property, we end up with the following equation for all k ≤ N:δg : ∂L ∂v (g k , v k+1 ) dτ k+1 = p T k dτ k Ad d k + h. ∂L ∂g (g k , v k+1 )(15)

Acknowledgements

The authors wish to thank the reviewers for their helpful feedback. This work was partly funded by the French ANR project SoHuSim and the French ANR project Morpho.


Jean-François Coulombel
email: jean-francois.coulombel@univ-nantes.fr

On the strong stability of finite difference schemes for hyperbolic systems in two space dimensions

Keywords: AMS subject classification: 65M12, 65M06, 35L45 Hyperbolic systems, finite difference schemes, stability

.

Introduction

Finite difference schemes are commonly used to approximate solutions to hyperbolic systems of conservation laws. In this article, we are interested in the stability of such finite difference schemes when applied to constant coefficients symmetric hyperbolic systems in two space dimensions. Symmetry is often crucial to show stability for a finite difference scheme, see for instance [START_REF] Yamaguti | Some remarks on the Lax-Wendroff finite-difference scheme for nonsymmetric hyperbolic systems[END_REF], and we shall therefore restrict to this framework.

For linear schemes, stability can be analyzed by means of Fourier transform and is found to be equivalent to uniform power boundedness of the so-called amplification matrix, see for instance [6, chapter 6]. The latter condition reads:

sup n∈N sup (ξ,η)∈R 2 C(ξ, η) n < +∞ ,
where we use from now on the notation (ξ, η) for the frequencies associated with the space variables (x, y) ∈ R 2 , and • denotes the spectral norm associated with the Hermitian norm for vectors:

∀ M ∈ M d (C) , M = sup X∈C d ,|X|=1 |M X| , |X| 2 := |X 1 | 2 + • • • + |X d | 2 .
One crucial ingredient in the analysis is the fact that the amplification matrix C(ξ, η) is a complex symmetric matrix, which simplifies the computation of its norm and/or its numerical radius. Two main subclasses of numerical schemes occur in practice:

1. Strongly stable schemes, for which the amplification matrix satisfies C(ξ, η) ≤ 1 for all (ξ, η) ∈ R 2 . This terminology dates back (at least) to [START_REF] Strang | On the construction and comparison of difference schemes[END_REF]. These are exactly the schemes for which the 2 -norm of the discrete solution decreases at each time step.

2. Schemes for which the numerical radius of C(ξ, η) is not larger than 1 for all (ξ, η). This property yields the bound sup

n∈N sup (ξ,η)∈R 2 C(ξ, η) n ≤ 2 .
The numerical radius is a very efficient tool for proving stability in specific situations, and this technique goes back to the seminal work [START_REF] Lax | Difference schemes for hyperbolic equations with high order of accuracy[END_REF]. The reader is also referred to the nice review [START_REF] Goldberg | On the numerical radius and its applications[END_REF] for a detailed introduction with further references. Strong stability is further studied in [START_REF] Tadmor | Complex symmetric matrices with strongly stable iterates[END_REF] where an equivalent condition is found in terms of real vectors. The latter condition is applied to the Lax-Wendroff scheme with stabilizer, and we shall refine some of these results below.

The main question we ask in this article is : does there exist a symmetric hyperbolic system with a numerical scheme that belongs to the second class without being strongly stable ? This seems far from obvious since as time goes by, one sometimes ends up showing that stability is in fact equivalent to strong stability. This is what we prove below for the so-called Lax-Wendroff scheme with stabilizer introduced in [START_REF] Lax | Difference schemes for hyperbolic equations with high order of accuracy[END_REF]. Even though strong stability may not be always equivalent to stability, it might still be recovered for sufficiently small CFL parameters. Quoting [12, page 128]: "It is interesting to speculate whether one can strengthen the basic conjecture by even claiming that when a symmetric hyperbolic scheme is stable, it is powerbounded with constant 2 and for sufficiently small ∆t/∆x even strongly stable." This is unfortunately not true, as detailed on a specific example below. To the best of our knowledge, this seems to be the first example of a numerical scheme of the second class that is shown not to be strongly stable. Still, our counterexample is not dissipative in the sense of Kreiss [START_REF] Kreiss | On difference approximations of the dissipative type for hyperbolic differential equations[END_REF]. In this more restrictive framework, we shall show that strong stability can be recovered for sufficiently small CFL parameters.

Let us conclude by observing that strong stability is also a powerful tool in the theory of discretized initial boundary value problems, see [START_REF] Coulombel | Semigroup stability of finite difference schemes for multidimensional hyperbolic initial-boundary value problems[END_REF]. We plan to investigate the influence of the absence of strong stability on semigroup estimates in a near future (this was the main motivation for studying the gap between the above two classes).

Main results

We consider a symmetric hyperbolic system in two space dimensions:

∂ t u + A ∂ x u + B ∂ y u = 0 , t ≥ 0 , (x, y) ∈ R 2 , u | t=0 = u 0 . (1)
The matrices A , B belong to M d (R) and are symmetric, so that the Cauchy problem (1) is well-posed in L 2 (R 2 ), see e.g. [START_REF] Benzoni-Gavage | Multidimensional hyperbolic partial differential equations[END_REF]. Moreover, the solution to (1) satisfies

∀ t ≥ 0 , u(t) L 2 (R 2 ) = u 0 L 2 (R 2 ) .
We introduce a finite difference approximation of (1). Let ∆x > 0 and ∆y > 0 denote some space steps in the x and y directions, and let ∆t denote the time step. Then the vector u n j,k , where (n, j, k) ∈ N × Z × Z, denotes an approximation of u(n ∆t, j ∆x, k ∆y). We define the CFL parameters

λ := ∆t ∆x , µ := ∆t ∆y ,
and for later use, we define the matrices

A := λ A , B := µ B .
We refer to [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]chapter IV.3], and [6, chapter 6] for a general description of finite difference schemes for two-dimensional hyperbolic systems, and we shall thus assume that the reader is familiar with the basic 2 -stability theory of finite difference schemes. In this article, we shall study the stability of three specific schemes:

• The upwind scheme:

u n+1 j,k = u n j,k - 1 2 A (u n j+1,k -u n j-1,k ) - 1 2 |A| (2 u n j,k -u n j+1,k -u n j-1,k ) - 1 2 B (u n j,k+1 -u n j,k-1 ) - 1 2 |B| (2 u n j,k -u n j,k+1 -u n j,k-1 ) . (2)
• The Lax-Wendroff scheme with stabilizer (see [START_REF] Lax | Difference schemes for hyperbolic equations with high order of accuracy[END_REF]):

u n+1 j,k = u n j,k - 1 2 A (u n j+1,k -u n j-1,k ) - 1 2 B (u n j,k+1 -u n j,k-1 ) - 1 2 A 2 (2 u n j,k -u n j+1,k -u n j-1,k ) - 1 2 B 2 (2 u n j,k -u n j,k+1 -u n j,k-1 ) + 1 8 (A B + B A) (u n j+1,k+1 -u n j+1,k-1 -u n j-1,k+1 + u n j-1,k-1 ) (3)
- 1 8 (A 2 + B 2 ) (4 u n j,k -2 u n j,k+1 -2 u n j,k-1 -2 u n j+1,k + u n j+1,k+1 + u n j+1,k-1 -2 u n j-1,k + u n j-1,k+1 + u n j-1,k-1
) .

• The Lax-Wendroff scheme without stabilizer (see again [START_REF] Lax | Difference schemes for hyperbolic equations with high order of accuracy[END_REF]):

u n+1 j,k = u n j,k - 1 2 A (u n j+1,k -u n j-1,k ) - 1 2 B (u n j,k+1 -u n j,k-1 ) - 1 2 A 2 (2 u n j,k -u n j+1,k -u n j-1,k ) - 1 2 B 2 (2 u n j,k -u n j,k+1 -u n j,k-1 ) (4)
+ 1 8 (A B + B A) (u n j+1,k+1 -u n j+1,k-1 -u n j-1,k+1 + u n j-1,k-1 ) .
We recall that in (2) the matrices |A|, |B| are defined as follows: let P, Q denote orthogonal matrices that diagonalize A and B:

P -1 A P = diag (α 1 , . . . , α d ) , Q -1 B Q = diag (β 1 , . . . , β d ) .
Then the matrices |A|, and |B|, are given by:

P -1 |A| P := diag (|α 1 |, . . . , |α d |) , Q -1 |B| Q = diag (|β 1 |, . . . , |β d |) .
Our main results are the following: Theorem 1.

• The scheme (2) is stable in 2 (Z 2 ) if and only if

∀ X ∈ R d , λ X, |A | X + µ X, |B | X ≤ |X| 2 ,
and in that case, it is even strongly stable.

• The scheme

(3) is stable in 2 (Z 2 ) if and only if ∀ X ∈ R d , λ 2 |A X| 2 + µ 2 |B X| 2 ≤ 1 2 |X| 2 ,
and in that case, it is even strongly stable.

The scheme (2) is a particular case of the finite volume scheme studied in [START_REF] Vila | Convergence of an explicit finite volume scheme for first order symmetric systems[END_REF]. It is obtained by choosing rectangles [j∆x 1 , (j + 1)∆x 1 ] × [k∆x 2 , (k + 1)∆x 2 ] as control volumes. For such control volumes, the stability condition obtained in [START_REF] Vila | Convergence of an explicit finite volume scheme for first order symmetric systems[END_REF]

λ + µ) max A , B ≤ 1 2 .
Our sufficient stability condition in Theorem 1 is less restrictive, but is is restricted to uniform cartesian grids while the condition in [START_REF] Vila | Convergence of an explicit finite volume scheme for first order symmetric systems[END_REF] applies to general triangulations.

The second part of Theorem 1 improves earlier criteria for the strong stability of (3). In [START_REF] Lax | Difference schemes for hyperbolic equations with high order of accuracy[END_REF], Lax and Wendroff prove stability for (3) under the assumption of Theorem 1 by means of the numerical radius. Abarbanel and Gottlieb [START_REF] Abarbanel | On proving stability of multidimensional difference schemes[END_REF] show strong stability under the assumption max( A , B ) ≤ 1/2 and, later on, Tadmor [START_REF] Tadmor | Complex symmetric matrices with strongly stable iterates[END_REF] shows strong stability under the assumption |A u| 4 +|B u| 4 ≤ 1/8 for all unit vector u. Theorem 1 gives a final optimal version to the stability analysis of (3).

In view of Theorem 1 and other known results in the literature, see e.g. [START_REF] Vaillancourt | On the stability of Friedrichs' scheme and the modified Lax-Wendroff scheme[END_REF][START_REF] Wendroff | Well-posed problems and stable difference operators[END_REF][START_REF] Zwas | On two step Lax-Wendroff methods in several dimensions[END_REF] for schemes in which the amplification matrix is either normal or a product of normal matrices, it might seem that strong stability is a common feature of many numerical schemes. The scheme (4) seems to be an exception for no strong stability result has been obtained so far in the case of noncommuting matrices. This is explained by our second main result. Theorem 2. Let A and B be given by A := 1 0 0 0 , B := 0 1 1 0 .

Then the scheme (4) is never strongly stable, though it is stable if λ 2/3 + µ 2/3 ≤ 1.

If A and B are both invertible, then the scheme (4) is strongly stable if the CFL parameters λ, µ are sufficiently small.

The paper is organized as follows. In section 3, we prove Theorem 1. Strong stability is proved directly on the amplification matrix for the scheme (2) while we use the decomposition into real and imaginary parts of the amplification matrix for the scheme (3). This gives us the opportunity to give a slight refinement of the results in [START_REF] Tadmor | Complex symmetric matrices with strongly stable iterates[END_REF]. Section 4 is devoted to the analysis of the scheme (4) in the particular case of the matrices given in Theorem 2. We shall also address the case of invertible matrices.

Strong stability results

The upwind scheme

The amplification matrix of the upwind scheme (2) can be written as

C(ξ, η) = I -2 sin 2 ξ 2 |A| + sin 2 η 2 |B| -i (sin ξ A + sin η B) .
The sufficient condition for stability is then obtained by choosing ξ = η = π. We now show that this condition is also sufficient for strong stability. We decompose the amplification matrix as follows:

C(ξ, η) = I -|A| -|B| + e iξ |A| -A 2 + e -iξ |A| + A 2 + e iη |B| -B 2 + e -iη |B| + B 2 .
All five real symmetric matrices

I -|A| -|B| , |A| -A 2 , |A| + A 2 , |B| -B 2 , |B| + B 2 ,
are nonnegative and their sum equals the identity matrix. The conclusion follows from the following general result1 whose proof is recalled for the sake of completeness:

Lemma 1. Let H 1 , . . . , H q ∈ M d (C) be nonnegative Hermitian matrices such that

H 1 + • • • + H q = I .
Then for all complex numbers z 1 , . . . , z q satisfying |z j | ≤ 1 for all j, there holds

z 1 H 1 + • • • + z q H q ≤ 1 .
Proof. We write

      z j H j 0 • • • 0 0 0 . . . . . . . . . 0 • • • • • • 0       = M * diag (z 1 I, . . . , z q I) M , M :=    H 1/2 1 0 • • • 0 . . . . . . . . . H 1/2 q 0 • • • 0    ,
where H

1/2 j denotes the unique nonnegative Hermitian square root of H j . The assumptions on the H j 's yield M = M * ≤ 1 and the result follows.

The spectral norm of complex symmetric matrices

The aim of this paragraph is to make a little more precise the statement of [START_REF] Tadmor | Complex symmetric matrices with strongly stable iterates[END_REF]Lemma 4.1] on the characterization of complex symmetric matrices with strongly stable iterates. Our result is:

Proposition 1. Let C ∈ M d (C
) be a symmetric matrix that we decompose as C = I -K + i J, where K and J are real symmetric matrices. Then C ≤ 1 if and only if for all real unit vectors x, y in R d there holds

K x, x K y, y + J x, y 2 ≤ 2 K x, x . (5)
Lemma 4.1 in [START_REF] Tadmor | Complex symmetric matrices with strongly stable iterates[END_REF] covers the "if" part of Proposition 1, and we briefly indicate why (5) gives a complete characterization of symmetric matrices in the unit ball of M d (C).

Proof. We start with the formula (3.1) in [START_REF] Tadmor | Complex symmetric matrices with strongly stable iterates[END_REF]. If C is a complex symmetric matrix whose real and imaginary parts are denoted R and J (that is, R = I -K with the notation of Proposition 1), there holds2 C = max

(u,v)∈R d ,|u| 2 +|v| 2 =1 R u, u -R v, v + 2 J u, v . (6)
Let x, y ∈ R d be unit vectors and let θ ∈ R. Then u := cos θ x and v := sin θ y satisfy

|u| 2 + |v| 2 = 1, so (6) gives 2 C ≥ 2 cos 2 θ R x, x -2 sin 2 θ R y, y + sin(2 θ) (2 J x, y ) = R x, x -R y, y + cos(2 θ) R x, x + R y, y + sin(2 θ) (2 J x, y ) .
Maximizing first over θ, and then over x, y, we get

2 C ≥ max (x,y)∈R d ,|x|=|y|=1 R x, x -R y, y + ( R x, x + R y, y ) 2 + 4 J x, y 2 .
The opposite inequality is shown by considering a couple (u, v) ∈ R d , |u| 2 + |v| 2 = 1, for which the maximum in ( 6) is attained. Such a couple can always be written under the form u := cos θ x, v := sin θ y, where x and y are unit vectors. Hence Corollary 3.2 in [START_REF] Tadmor | Complex symmetric matrices with strongly stable iterates[END_REF] can be improved as

2 C = max (x,y)∈R d ,|x|=|y|=1 R x, x -R y, y + ( R x, x + R y, y ) 2 + 4 J x, y 2 . (7)
Recalling [START_REF] Tadmor | Complex symmetric matrices with strongly stable iterates[END_REF]Lemma 4.1], it remains to prove that if C is a complex symmetric matrix whose norm does not exceed 1, then (5) holds. Let us first observe that for all real unit vector x, the complex number C x, x belongs to the unit disk, so its real part R x, x belongs to [-1, 1]. In other words, there holds R ≤ 1. Since C ≤ 1, we can apply [START_REF] Kreiss | On difference approximations of the dissipative type for hyperbolic differential equations[END_REF] and obtain ( R x, x + R y, y ) 2 + 4 J x, y 2 ≤ 2 -R x, x + R y, y , for all unit vectors x, y (the right hand side is nonnegative since R ≤ 1). Squaring the inequality and simplifying some terms, we obtain

R x, x 2 + J x, y 2 ≤ 1 -( R x, x -R y, y ) (1 -R x, x ) ,
which is nothing but (5) by recalling R = I -K.

The Lax-Wendroff scheme with stabilizer

Our goal is now to prove that the Lax-Wendroff scheme with stabilizer (3) is stable, and even strongly stable, if and only if 2 (A 2 + B 2 ) ≤ I. The amplification matrix of (3) reads C(ξ, η) = I -K + i J with

J := sin ξ A + sin η B , K := 1 2 J 2 + 2 sin 2 ξ 2 + sin 2 η 2 sin 2 ξ 2 A 2 + sin 2 η 2 B 2 . (8)
Following [START_REF] Lax | Difference schemes for hyperbolic equations with high order of accuracy[END_REF], we choose ξ = η = π and find that the condition 2 (A 2 + B 2 ) ≤ I is necessary for stability. Let us now show that it is also sufficient for proving strong stability of (3). In view of Proposition 1, see also [START_REF] Tadmor | Complex symmetric matrices with strongly stable iterates[END_REF]Corollary 4.3], it is sufficient to show that for all real unit vectors x, y in R d there holds

3 K x, x K y, y ≤ (2 K -J 2 ) x, x . (9)
We use the definition [START_REF] Lax | Difference schemes for hyperbolic equations with high order of accuracy[END_REF] and thus wish to show the inequality

K x, x K y, y ≤ 4 sin 2 ξ 2 + sin 2 η 2 sin 2 ξ 2 |A x| 2 + sin 2 η 2 |B x| 2 . ( 10
)
For simplicity, we use the notation

α := sin 2 ξ 2 , β := sin 2 η 2 .
We estimate K y, y for all real unit vector y. The matrix K reads

K = 2 α A 2 + 2 β B 2 + 2 α β (A 2 + B 2 ) + 1 2 sin ξ sin η (A B + B A) . Since |A y| 2 + |B y| 2 ≤ 1/2, we can write |A y| 2 = r 2 (1 + sin θ) , |B y| 2 = r 2 (1 -sin θ) , 0 ≤ r ≤ 1 2 , |θ| ≤ π 2 .
We compute

K y, y = r (α + β) + 2 r α β + r sin θ (α -β) + sin ξ sin η A y, B y ≤ r (α + β + 2 α β) + r sin θ (α -β) + | sin ξ| | sin η| |A y| |B y| ,
and we thus derive the estimate

K y, y ≤ r (α + β + 2 α β) + r sin θ (α -β) + cos θ | sin ξ| | sin η| 2 ≤ r (α + β + 2 α β) + r (α -β) 2 + sin 2 ξ sin 2 η 4 1/2 . It remains to compute (α -β) 2 + sin 2 ξ sin 2 η 4 1/2 = (α + β -2 α β) 2 1/2 = α + β -2 α β ,
where the last equality comes from the fact that α, β ∈ [0, 1]. Eventually we have derived the estimate K y, y ≤ 2 r (α + β) ≤ α + β .

Coming back to [START_REF] Strang | On the construction and comparison of difference schemes[END_REF], we see that it only remains to show

K x, x ≤ 4 α |A x| 2 + β |B x| 2 ,
which is precisely what is obtained in [11, page 75]. Let us recall the derivation of this final estimate for the sake of clarity. The definition (8) gives

K x, x = 1 2 |J x| 2 + 2 (α + β) α |A x| 2 + β |B x| 2 .
The norm of the vector J x is estimated as follows: We have thus shown that (10) holds, and Proposition 1 shows that the scheme (3) is strongly stable.

might be recovered for sufficiently small CFL parameters. This assertion should be taken very cautiously though, and we have no other example of this fact except (4).

Eventually, let us observe that strong stability can hold for (4) with noncommuting matrices, even when the CFL parameters are not small. Let us consider for instance the case A := 1 0 0 -1 , B := 0 1 1 0 , which corresponds to the two-dimensional wave equation. Then the amplification matrix of (4) reads 1 -2 λ 2 sin 2 ξ 2 -2 µ 2 sin 2 η 2 I -i (sin ξ λ A + sin η µ B ) .

This is a normal matrix whose spectral radius is not larger than 1 for all (ξ, η) if and only if λ 2 +µ 2 ≤ 1. In that case, stability is equivalent to strong stability. As expected from the general theory in [START_REF] Turkel | Symmetric hyperbolic difference schemes and matrix problems[END_REF], the stability domain encompasses the "scalar" one (λ 2/3 + µ 2/3 ≤ 1).

|J x| 2 ≤

 2 sin 2 ξ |A x| 2 + sin 2 η |B x| 2 + 2 | sin ξ| | sin η| |A x| |B x| = 4 α (1 -α) |A x| 2 + 4 β (1 -β) |B x| 2 + 8 α (1 -β) β (1 -α) |A x| |B x| ≤ 4 α (2 -α -β) |A x| 2 + 4 β (2 -α -β) |B x| 2 ,which eventually gives the expected estimate:K x, x ≤ 4 α |A x| 2 + β |B x| 2 .

This result seems quite classical. It appears as Exercise 301 in the additional list of exercises of[START_REF] Serre | Matrices. Graduate Texts in Mathematics[END_REF], see http://www.umpa.ens-lyon.fr/∼serre/DPF/exobis.pdf. It may probably be found in earlier textbooks or references that the author is not aware of.

[START_REF] Benzoni-Gavage | Multidimensional hyperbolic partial differential equations[END_REF] The equality can also be written, as in[START_REF] Tadmor | Complex symmetric matrices with strongly stable iterates[END_REF], with 2 | J u, v | on the right hand side by changing u into -u, but this is of no consequence.

This inequality implies (5) thanks to Cauchy-Schwarz.

Acknowledgements Research of the author was supported by the French Agence Nationale de la Recherche, contract ANR-08-JCJC-0132-01.

The Lax-Wendroff scheme without stabilizer

Let us first recall that the scheme ( 4) is known to be stable under the condition that for all real unit vector u, there holds |A u| 2/3 + |B u| 2/3 ≤ 1, see [START_REF] Turkel | Symmetric hyperbolic difference schemes and matrix problems[END_REF]. (The reader can also consult [START_REF] Goldberg | On the numerical radius and its applications[END_REF] for the more restrictive -though easier -criterion |A u| 2 + |B u| 2 ≤ 1/4.)

We first study the finite difference scheme (4) when the matrices A and B are given by:

The amplification matrix of the Lax-Wendroff scheme (4) reads C(ξ, η) = I -K + i J with

Let us assume that the scheme is strongly stable. In particular, for all η ∈ R, there holds C(π, η) ≤ 1. Choosing x := e 2 and y := e 1 the vectors that span the canonical basis of R 2 , Proposition 1 gives e 2 , K e 2 e 1 , K e 1 + J e 2 , e 1 2 ≤ 2 e 2 , K e 2 , that is,

The latter condition reads

which gives obviously a contradiction by choosing η > 0 small enough.

Let us now assume that the matrices A and B are invertible, and let us prove that for sufficiently small CFL parameters λ, µ the amplification matrix defined by [START_REF] Tadmor | Complex symmetric matrices with strongly stable iterates[END_REF] has strongly stable iterates. Our goal is to show that the criterion [START_REF] Serre | Matrices. Graduate Texts in Mathematics[END_REF], which is sufficient for strong stability, is satisfied for λ, µ small enough. We still use the notation

From the definitions (11), we compute

where we have used the invertibility of A and B , and the constant c > 0 is independent of α, β, λ, µ and x. Let us now give a bound for K :

where C denotes a positive constant, which may vary from one line to the next, that does not depend on α, β, λ, µ. Hölder's inequality yields

so choosing λ, µ small enough, the inequality ( 9) is satisfied. This tends to indicate that if one restricts to numerical schemes that are dissipative in Kreiss's sense, then strong stability