Blandine Bolteau

Kevin Dos Santos

Frédéric Gelebart

José Gomes

Jérémie Teisseire

Etienne Barthel

Jérôme Fresnais
email: jerome.fresnais@sorbonne-universite.fr

A new platform for gravitational microfluidic using ferrofluids

Keywords: superhydrophobic surfaces, ferrofluids, microfluidic, tilted surface

Among the large variety of microfluidic platforms, surface devices are a world apart. Electrowetting systems are used to control the displacement of droplets among predetermined pathways. More confidential, superhydrophobic surfaces are more and more described as new elements to guide spherical droplet reactors. As such, they can exhibit confinement properties analog to channel-based microfluidics.

In this article, we describe a new strategy to use superhydrophobic surfaces as permanently-tilted microfluidic platform, on which droplets containing iron oxide nanoparticles are guided with permanent magnets. These droplets are fed with water through a capillary tube until their weight exceeds the magnetic field force. Thus, the volume at which the droplet rolls off the surface is only governed by the initial quantity of magnetic nanoparticles and the tilting angle of the surface. This phenomenon provides a strategy for droplet dilution in a simple and reproducible manner, which is not that easy in micro-channels, and a key advantage of open systems.

As a proof of concept, we used this platform to prepare magnetic filaments by a salting-out process already described in large batches. By reducing salt concentration on the platform, we are able to control the electrostatic attractive interactions between iron oxide nanoparticles coated with poly(acrylic acid) and a positively charged polyelectrolyte (poly(diallyldimethylammonium chloride)). The formation of nanostructured filaments was conducted in two minutes while more than 30 minutes were required by dialysis. Our results also illustrate the power of microfluidic reaction processes since such magnetic filaments could not be obtained through direct batch dilution due to mixing issues. Such microfluidic platforms could be useful for the efficient and simple dilution of systems where reactivity is controlled by concentration.

Introduction

In the framework of microfluidic, open-channel systems are less used than closed microchannels, even-though their advantages are numerous. They exhibit the potential of confined liquid systems without the detrimental interferences with the micro-channel walls. In this area, electro-wetting systems are interesting because electric field application allows the guidance of droplets on a predetermined pathway [START_REF] Pollack | Electrowetting-Based Actuation of Droplets for Integrated Microfluidics[END_REF] . Platforms are described to use trains of droplets that can merge or separate on-demand. More confidential is the use of superhydrophobic surfaces as microfluidic open-channel systems on their own [START_REF] Seo | A Superhydrophobic Magnetic Elastomer Actuator for Droplet Motion Control[END_REF][START_REF] Oliveira | Two-Dimensional Open Microfluidic Devices by Tuning the Wettability on Patterned Superhydrophobic Polymeric Surface[END_REF] . There are some examples where droplets are guided when elastic surfaces are deformed [START_REF] Nagappan | Emerging Trends in Superhydrophobic Surface Based Magnetic Materials: Fabrications and Their Potential Applications[END_REF] . Fewer examples are given for « ferrofluidics » [START_REF] Bormashenko | New Investigations on Ferrofluidics: Ferrofluidic Marbles and Magnetic-Field-Driven Drops on Superhydrophobic Surfaces[END_REF] . The advantages of this technique is that no specific development is required on the hydrophobic corrugated surface itself, and the driving force is low because of the reduced hysteresis of the superhydrophobic state [START_REF] Cassie | Wettability of Porous Surfaces[END_REF] . These ferrofluid droplets can be displaced and deformed by controlling the magnetic field amplitude and gradient [START_REF] Zhu | Nonlinear Deformation of a Ferrofluid Droplet in a Uniform Magnetic Field[END_REF] . Additionally, recent articles present the use of superhydrophobic magnetic needle to handle magnetic droplets [START_REF] Yang | Automatic Magnetic Manipulation of Droplets on an Open Surface Using a Superhydrophobic Electromagnet Needle[END_REF] . Here we propose a new route to use superhydrophobic surfaces as platforms for of milli/microfluidic using ferrofluids not only as mobile droplets but also as reactive ones. Ferrofluid droplets can be moved on surfaces using permanent magnets. These droplets can be fed with solutions containing reactive entities (either molecules, ligands, polymers, proteins, or even cells) that can react (or not) with the magnetic nanoparticles. The dilution that occurs during the droplet feeding is well controlled and is a valuable manner to progressively reach the targeted concentration inside the droplets, which is not easily achieved in closed channels for large dilution rates, except on electro-wetting droplet platforms. Moreover, a very simple way to obtain reproducible mixing ratios between the on-growing droplet and the reactive feeding is to tilt the surface so that the increasing weight of the droplet eventually reaches a threshold where it overcomes the magnetic confinement. Evidences have been demonstrated in the literature where tilting the surface can move the droplet from one hydrophilic post to another one [START_REF] Kong | Motorized Actuation System to Perform Droplet Operations on Printed Plastic Sheets[END_REF] . The balance between both forces (weight and magnetic coupling) is only driven by the ferrofluid concentration, the starting droplet volume, and the tilting angle of the surface. Thus, reproducible dilution / mixing ratios can be achieved and potentially parallelized using networks of magnets.

Methods

Superhydrophobic surfaces

Superhydrophobic surfaces were prepared by a soft lithography process using a SU-8 photosensitive resin mold. It consists of a square network of cylindrical pillars with a diameter of 17.1 µm, 61.6 µm high and a spacing of 50 µm between the center of pillars. The SU-8 resin was deposited through spin coating (30 s at 500 rpm, followed by 30 s at 1900 rpm) and preheated at 65°C during 3 minutes and 95°C during 9 minutes. Exposure through a transparent mask was carried out with 8 cycles of 30 seconds on and 20 seconds off at 70 % of the maximum power of the lamp at 365 nm in a UV-KUBE. The resin was post-baked at 65°C during 2 minutes followed by one 7 minutes at 95°C to ensure a complete polymerization of the resin. Development is conducted with SU8-développeur (Propylene Glycol Monomethyl Ether Acetate, PG-MEA). At last, a hard-bake is conducted at 150°C during 15 minutes. Then, SYLGARD 184 was used to prepare a negative mold of the pillars. After activation with an O2 plasma treatment and a silanization with trichloro(1H,1H,2H,2Hperfluorooctyl)silane in the vapor phase, the mold was used to reproduce the pillar network using SYLGARD 184. It can be re-used several times for further replication of the structures. The surfaces were characterized by Scanning Electron Microscopy (SEM) and optical microscopy (OLYMPUS…). Contact angles (advancing and receding), as well as sliding angles were measured with a portable goniometer (Krüss DSA 100) and a home-made tilting goniometer. Distilled water was used with droplet volumes ranging from 5 µL to 20 µL. ImageJ software with Snake plugin was used to measure the contact angles of droplets deposited on the surfaces.

Ferrofluid

The chemicals were used as received: Nitric acid (HNO3, 69.5% water solution, Carlo Erba), hydrochloric acid (HCl, 37% water solution, AnalaR Normapur), iron (II) chloride (FeCl2, 4H2O, AnalaR Normapur, VWR), iron (III) chloride (FeCl3, 6H2O, Prolabo), iron (III) nitrate (Fe(NO3)3, 9H2O, Technical, VWR), ammonia aqueous solution (NH4OH, 28-30%, Sigma-Aldrich). Water was high purity grade, obtained by an Elga purification system (18 MΩ cm). Aqueous solutions of maghemite MNP (γ-Fe2O3) were synthesized following Massart's coprecipitation method [START_REF] Bee | Synthesis of Very Fine Maghemite Particles[END_REF] . Ammonia was added to an equimolar acidic solution of iron (II) and iron (III) to induce the co-precipitation into Fe 3 O 4 (magnetite). It was then boiled for 30 minutes in the presence of iron nitrate in order to oxidize the magnetite into maghemite (γ-Fe 2 O 3 ). The resulting dispersion was then washed with acetone and ether before redispersion in an aqueous nitric acid solution.

Iron oxide nanoparticles were coated with poly(acrylic acid) (PAA) using a well-established process [START_REF] Sehgal | Precipitation-Redispersion of Cerium Oxide Nanoparticles with Poly(Acrylic Acid): Toward Stable Dispersions[END_REF][START_REF] Guibert | Controlling Nanoparticles Dispersion in Ionic Liquids by Tuning the PH[END_REF] . In brief, a dispersion of iron oxide at 0.5 %wt. (100 mL) is added to an acidic solution of PAA at 1 % wt. (300 mL) and stirred 20 minutes with the use of a magnetic bar. The particles spontaneously aggregate due to their interactions with uncharged PAA. These aggregates are magnetically separated and the supernatant is removed. The pH of the remaining preparation is raised to 10. At this pH value, the nanoparticles are well dispersed and stable in a pH range of 10 to 5.

Those nanoparticles can be aggregated through their interactions with positive polyelectrolytes. Poly(dialyldimethylammonium chloride) (PDADMAC) was used in this study [START_REF] Yan | Magnetic Nanowires Generated via the Waterborne Desalting Transition Pathway[END_REF] . To control the attractive electrostatic interactions between nanoparticles and positive polyelectrolyte, brine solution of PDADMAC in 0.5 M NH 4 Cl and a dispersion of PAA-coated nanoparticles at the same ionic strength were prepared. Mixing the two salty solutions, a stable dispersion is obtained for which hydrodynamic size reflects well dispersed PAA-coated nanoparticles. Decreasing the ionic strength through dilution down to 0.1 M, aggregates form [START_REF] Fresnais | Nanoparticle Aggregation Controlled by Desalting Kinetics[END_REF] . If this decrease in conducted in the presence of a magnetic field, the aggregates are elongated [START_REF] Fresnais | Electrostatic Co-Assembly of Iron Oxide Nanoparticles and Polymers: Towards the Generation of Highly Persistent Superparamagnetic Nanorods[END_REF] .

Ferrofluid droplet guidance

To move ferrofluid droplets on a tilted superhydrophobic surface, magnets was attached onto an oscillating motor that can move with a rotation angle of 180° max and adjustable tangential velocity ranging from 0.01 cm/s up to 4 cm/s. The motor is driven by an Arduino type set-up and programmed with the Arduino IDE software. The magnets consist of a column of 2 cylindrical magnets with 1 cm in diameter and 2 cm in height and a conical magnet on top (large diameter is 1cm, small diameter is 0.5 cm) made of Neodyn-Bore magnet (SUPERMAGNETE). The distance between magnet and surface (i.e. bottom of the droplet) varies between 2 and 20 mm.

Results & Discussion

Measurements with water droplets At As a first step, PDMS superhydrophobic surfaces (image and structure characterization are given in figure 1) were synthesized according to classical soft photolithography processes [START_REF] Du Roure | Microfabricated Arrays of Elastomeric Posts to Study Cellular Mechanics[END_REF] . These surfaces are made of arrays of pillars ca. 17 µm in diameter and 63 µm high. They were characterized in terms of contact angle, contact angle hysteresis, and sliding angle using water. The advancing contact angle was measured at 165°, the receding contact angle was measured at 145° by tilting the surface until the droplet moves. These two contact angles give an hysteresis of 20°. The sliding angles were measured with water droplets from 5 µL to 20 µL. As expected, the smaller the droplets, the higher the sliding angle of the surface (figure 2). These values are in good agreement with the results from the literature on those such model surfaces [START_REF] Joanny | A Model for Contact Angle Hysteresis[END_REF][START_REF] Dubov | Superhydrophobic Silica Surfaces: Fabrication and Stability[END_REF] . Indeed, the sliding angles can be fitted comparing the hysteresis value and the weight of the droplet as

mg sin α > Rγ(cos θ r -cos θ a ), ( 1
)
where m is the weight of the droplet, g is the gravitational acceleration, R is the diameter of the droplet, γ is the surface tension of the liquid and θ r and θ a are the receding and advancing contact angles. We found no significant difference for the sliding angles with the ferrofluid droplets which is reasonable as surface tension changes are minute between the different dispersions and water. These measurements are important to evaluate the role of viscous friction for droplet motion. A simple calculation based on Eq. 1, for a 5 µL droplet, gives a friction force of the order of 10 -2 µN, which is negligible compared to the magnetic force applied to the ferrofluid droplet by the magnet (see calculation below).

Measurement with ferrofluid droplets

Without magnetic forces applied to the droplet, Eq. 1 makes it straightforward to evaluate the tilting angle at which a droplet will roll off. On the contrary, when a magnet is used to confine a ferrofluid droplet, the tilting angle at which the magnetic droplet can be manipulated is larger than the roll-off angle without magnet. Typically, tilting angles of 10°, 20°, 40°, and 60° will be used with ferrofluid droplet volumes ranging from 5 to 20 µL.

To evaluate the ideal distance between the magnet and the surface to move ferrofluid droplets, we investigated the form of the droplet (figure 3) versus the distance between the top of the magnet and the surface. The deformation of the droplet has already been established on ferrofluid droplets with magnetic fields comparable to ours. Indeed, the maximum magnetic field applied in a previous study [START_REF] Zhu | Nonlinear Deformation of a Ferrofluid Droplet in a Uniform Magnetic Field[END_REF] is 370 mT, which is in good agreement with our measurements which give a magnetic field amplitude of 335 mT at 2 mm from the magnet. For such a magnetic field and magnetic field gradient amplitudes, the magnetic force applied to the ferrofluid droplet is:

(

where corresponds to the difference of magnetic susceptibility of the ferrofluid, m is the mass of magnetic materials, is the density of the fluid, , B and gradB are the magnetic field and magnetic gradient amplitude respectively. In our case, this magnetic force is of the order of 1 µN, which is 100 times larger than the typical droplet friction forces. When the magnet approaches the droplet, the droplet flattens as observed by Zhu et al [START_REF] Zhu | Nonlinear Deformation of a Ferrofluid Droplet in a Uniform Magnetic Field[END_REF] . Additionally, we also measure the contact angle as the magnet gets closer to the ferrofluid droplet. Indeed, for the droplet to be driven by the magnet, the transition to the Wenzel state must be avoided. Figure 4 represents the contact angle of ferrofluid droplets with different concentrations ranging from 0.4 %wt. to 6.4 %wt. We find that the contact angle starts around the mean contact angle (ca. 145°) for large magnet-droplet distances. As the magnet approaches and the droplet flattens, the contact angle increases towards the advancing contact angle (ca. 165°-170°). In figure 4a, we have reported the variation of the relative aspect ratio of the droplet versus the distance from the magnet for various concentrations. The data are normalized by the aspect ratio of the droplet when the magnet is 20 mm away. The deformation is larger and appears at a larger distance as the concentration increases. As a consequence, the contact angle varies as the droplet is deformed, following the same dependence with ferrofluid concentration and surface to magnet distance. In contrast to literature reports, no transition toward the Wenzel state is observed in our case, even when the magnet is as close as 2 mm from the droplet. The good stability of the Cassie Baxter regime is not only due to the low surface fraction (9.4%) but also to the length of the pillars. Indeed, Al-Azawi et al. demonstrated the transition towards Wenzel state for surface fraction values 5 times larger and for pillars which are 40 µm high compared to 63 µm in the present study [START_REF] Al-Azawi | Friction and Wetting Transitions of Magnetic Droplets on Micropillared Superhydrophobic Surfaces[END_REF] . They reported that the transition towards the Wenzel state was enhanced by oscillations of the droplet on the surface, which were not observed in our experiments. In order to apply this system to a microfluidic platform, we evaluated the distance at which droplets can be manipulated with the magnet as a function of nanoparticle concentration (figure 5). The distance at which droplets follow the motion of the magnet depends linearly on the starting volume and the tilting angle for droplet concentration larger than 0.8%wt. Below this concentration, it is not possible to move the droplets, whatever the distance between the magnet and the droplet base. We thus chose to work with a ferrofluid concentration of 1%wt. and a distance between the magnet and the droplet of 4 mm, to ensure active and responsive droplets when the magnet is moved below the surface. We investigated the roll-off volume of a ferrofluid droplet moving periodically on a tilted surface. At each period, the droplet was fed with water through a capillary tube situated at the top of the droplet trajectory until the growing droplet escaped from the magnetic confinement and rolled off the surface. Results are reported in figure 6a (error bars correspond to the standard deviation of at least 3 measurements). An important observation is that the volume at roll-off normalized by the initial volume of the droplet is nearly constant. Thus, dilution at roll-off is independent from initial volume. Figure 6b reports dilution ratio versus the tilting angle for a ferrofluid with initial concentration of 1 %wt., extracted from figure 6a. Therefore, for a given initial ferrofluid concentration, dilution at roll-off strongly depends on tilting angle. With such a device, it is straightforward to select the dilution at droplet roll-off. To evaluate the potential of this method for rapid and robust droplet dilution, we investigated the volume at roll-off for a tilting angle of 40° and different droplet velocities ranging from 8 mm/s up to 37.5 mm/s (figure 7). The roll-off volumes are found to be superimposed for the different initial ferrofluid droplets, whatever the droplet velocity. These velocities can be compared to droplet generation and flows in micro-channel in biphasic systems where flow rates as high as 150 mm/s are measured. Thus, manipulation of droplets with our device can be considered for the dilution of numerous droplets in short times.

Figure 7: Total volume at roll-off versus initial droplet volume at different droplet velocities (8, 20, and 37.5 mm/s) for a tilting angle of 40° (at the concentration of 1 %wt.). The line corresponds to a linear fit.

Measurements with reactive magnetic droplets

As a proof of concept, we investigated the reaction of iron oxide nanoparticles towards electrostatic complexation with oppositely charged polyelectrolytes. Indeed, we have previously demonstrated a controlled pathway to tune the interaction between positively charged polyelectrolytes and negatively charged nanoparticles [START_REF] Fresnais | Nanoparticle Aggregation Controlled by Desalting Kinetics[END_REF] . Addition of a large quantity of salt (ca. NH 4 Cl at 0.5 M) in a dispersion of poly(acrylic acid) coated nanoparticles and in a positively-charged polyelectrolyte solution screens their interaction when they are mixed. Dilution with water decreases salt concentration below a threshold at which attractive electrostatic interactions are revealed. We have previously demonstrated that rods can be obtained in 1 hour through the dialysis of the salted mixture in the presence of a magnetic field [START_REF] Yan | Magnetic Nanowires Generated via the Waterborne Desalting Transition Pathway[END_REF] . Here we investigated the dilution process in a droplet of given initial volume guided by the magnet and diluted with water through a capillary tube. Figure 8b represents the picture of the droplet containing the magnetic nanoparticles as well as the positive polyelectrolyte and salt. Dilution occurs at every passage of the droplet above the capillary tube. Interestingly, the magnetic nanoparticles are not uniformly distributed in the droplet while water is added to the moving droplet. This is due to the strong magnetic interaction between the nanoparticles and the magnet. On the contrary, the non-magnetic species (salts, polyelectrolytes) can diffuse freely inside the droplet, which allows the formation of the aggregates at the interface where the salt concentration reaches the critical value (0.5 M of NH 4 Cl). This phenomenon will be studied in a forthcoming work. At the end of the dilution (the volume of the droplet is increased 10 times from its initial volume), the droplets are collected and their content is analyzed by optical (figure 8-c) and transmission electron microscopy (figure 8-d). As a comparison, image 8a represents the pristine nanoparticles to be co-assembled. We evidenced the synthesis of elongated aggregates as found in our previous work [START_REF] Yan | Magnetic Nanowires Generated via the Waterborne Desalting Transition Pathway[END_REF] . In previous attempts, these aggregates could not be obtained through a dilution process, due to the hydrodynamic perturbations that usually occur when two liquids are mixed together and also because of the rapid assembly rate between oppositely charged species above the critical salt concentration. This is thus the first evidence of the synthesis of those elongated aggregates through a dilution process. Moreover, assembly is obtained within minutes instead of a one-hour dialysis. Table 1 reports the mean lengths and their corresponding standard deviations for different starting droplet volumes from 5 to 30 µL. The bigger the droplet, the longer the rods. The differences observed in length are not that large but it could be possible to change drastically their length by changing the strength of the magnetic field (the distance between the droplet and the magnet) or the concentration of ferrofluid.

Conclusion

To conclude, we propose a simple and robust method to perform dilution in microfluidic systems using ferrofluid droplets moving on a tilted superhydrophobic surface. Control over progressive dilution is not easy to achieve within microfluidic systems. Here, dilution at rolloff only depends on the tilting angle and not on the initial droplet volume, and is very reproducible for various velocities up to 37.5 mm/s. This dilution process was successfully used to prepare nanostructured rods of ca. 10 µm long in a few minutes. Similar rods cannot be obtained by direct dilution of the bulk solution and their synthesis takes one hour through dialysis. This dilution set-up could be useful to achieve dilution of solutions where a solute can react with the surface of iron oxide nanoparticles. For example, investigation of kinetic adsorption of molecules of interest such as proteins or cancer markers in blood drops could be achieved rapidly and in a reproducible manner.

Figure 1 :

 1 Figure 1: Schematics with pillar geometry; SEM image.

Figure 2 :

 2 Figure 2: Sliding angle versus droplet volume; line corresponds to equation 1.

Figure 3 :

 3 Figure 3: Images of a ferrofluid droplet (0.8 M) at three decreasing magnet / surface distances (20 mm, 10 mm, and 2 mm).

  b

Figure 4 :

 4 Figure 4: (a) contact angle variation versus distance between the magnet and the droplet for different ferrofluid concentration (from 6.4 %wt. to 0.4 %wt.); (b) relative aspect ratio of the droplet versus the distance between the magnet and the droplet

Figure 5 :

 5 Figure 5: minimal distance between the magnet and the droplet versus nanoparticle concentration (%wt.) to move the droplet on the surface. The grey area corresponds to the

Figure 6 :

 6 Figure 6: (a) Ratio between roll-off volume and initial volume versus the initial droplet volume at different tilting angles (10°, 20°, 40°, and 60°) for a ferrofluid at 1 %wt.; (b) dilution factor versus tilting angle

Figure 8 :

 8 Figure 8: Image of a droplet containing a salted ferrofluid with positively charged polyelectrolytes (b). The capillary tube delivers water so that interparticle electrostatic interactions are revealed by salt dilution. Optical microscopy (c) and TEM images of the

Table 1 :

 1 magnetic rod mean lengths and standard deviation of rods versus initial droplet volume

	Initial droplet volume (µL)	Rods mean length (µm)	Standard deviation
	5	9.8	3.8
	10	11.8	3.6
	20	11	5.3
	30	13.2	6.3

Acknowledgements

[Insert text here]


I Charpentier
email: icharpentier@unistra.fr

Bruno Cochelin

Towards a full higher-order AD-based continuation and bifurcation framework

Keywords: nonlinear parametric problems, continuation method, bifurcation analysis, geometric series, higher order automatic differentiation, Diamant AMS Subject Classification: 41A58, 65P30, 68N19, 70E60

Some of the theoretical aspects of continuation and bifurcation methods devoted to the solution for nonlinear parametric systems are presented in a higher-order automatic differentiation (HOAD) framework. Besides benefits in terms of generality and ease of use, HOAD is used to assess fold and simple bifurcations points. In particular, the formation of a geometric series in successive Taylor coefficients allows for the implementation of an efficient detection and branch switching method at simple bifurcation points.

Some comparisons with the Auto and MatCont continuation software are proposed. Strengths are then exemplified on a classical case study in structural mechanics.

Introduction

Numerical continuation and bifurcation analysis are classical tools for the study of nonlinear parametric equations, see [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] and the references therein for instance. In general terms, these methods allow for the solution of nonlinear algebraic systems of n equations,

R(U ) = R(u, λ) = 0, (1)
where the unknown U = (u, λ) of the nonlinear function R comprises the state vector u of dimension n and some scalar control parameter λ to be varied. Frequently, this parameter has a physical meaning, representing a load parameter in structural mechanics or a volumetric flow rate in fluid mechanics, for instance. The solutions of the underdetermined problem (1) are one-dimensional continua, referred to as solution branches, that may contain singular points such as folds (limit points), and bifurcation points where two branches intersect. Note that a number of additional bifurcations may arise in the solution of nonlinear parametric ordinary differential equations [START_REF] Dhooge | New features of the software matcont for bifurcation analysis of dynamical systems[END_REF][START_REF] Eusebius | AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Eq February[END_REF][START_REF] Seydel | Numerical computation of branch points in nonlinear equation[END_REF]. Such singular points indicate a qualitative change in the behavior of the system under study. They can be detected and classified from determinant calculations and/or derivative computations, see for instance [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory -Volume I[END_REF][START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF].

During the last thirty years, general-purpose software, free and commercial ones, have been proposed to engineers and scientists to draw bifurcation diagrams without embarking into the heavy task of programming their own continuation algorithm. Most of these continuation software, see MatCont [START_REF] Dhooge | New features of the software matcont for bifurcation analysis of dynamical systems[END_REF] and Auto [START_REF] Eusebius | AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Eq February[END_REF] for instance, follow the first order predictor-corrector principles described in [START_REF] Keller | Lectures on numerical methods in bifurcation problems[END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]. Besides, under analyticity assumptions, solution branches may be approximated as higher-order truncated Taylor expansions [START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF][START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF][START_REF] Cochelin | Méthode asymptotique numérique[END_REF] to obtain continuous solutions and representations. Their range of validity is estimated a posteriori from the remainder of the series to provide large adaptive stepsize.

Today, Taylor-based nonlinear solvers take full benefit of higher-order automatic differentiation (HOAD), notably in terms of generality, efficiency and automation, see [START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF] and references therein. Our HOAD continuation framework, namely Diamant, moreover agrees with homotopy for the solution of nonlinear eigenvalue problems and with sensitivity analysis for both continuation and homotopy methods [START_REF] Charpentier | Sensitivity computations in higher order continuation methods[END_REF]. It has been recently adapted to the determination of singular boundaries [START_REF] Hentz | Higher-order continuation for the determination of robot workspace boundaries[END_REF] of particular piecewise smooth nonlinear systems by using the extended system of equations introduced in [START_REF] Seydel | Numerical computation of branch points in nonlinear equation[END_REF].

Higher order information can be used to locate singular points [START_REF] Van Dyke | Analysis and improvement of perturbation series[END_REF][START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory -Volume I[END_REF][START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF] and AD can be used to extract it [START_REF] Dhooge | New features of the software matcont for bifurcation analysis of dynamical systems[END_REF][START_REF] Guckenheimer | Computing periodic orbits and their bifurcations with automatic differentiation[END_REF]. In particular, simple folds can be detected in a series in straightforward manner. As discussed for higher-order perturbation methods [START_REF] Van Dyke | Analysis and improvement of perturbation series[END_REF], the formal power series arising near to a singularity reveals the analytical structure of the solution. It can be used to "improve the utility of the series" and the convergence of the numerical method. This has been applied to specific quadratic formulations of nonlinear parametric problems near to simple bifurcation points [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF]. In this paper, we take advantage of HOAD to rework this proposal in the general framework of (1) in order to implement an efficient detection and branch switching method.

Two tutorial examples from MatCont [START_REF] Kuznetsov | Tutorial II: One-parameter bifurcation analysis of equilibria with matcont[END_REF] and Auto [START_REF] Doedel | Numerical analysis and control of bifurcation problems (i) bifurcation in finite dimensions[END_REF][START_REF] Eusebius | AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Eq February[END_REF] are used to illustrate some strengths and limitations of the Matlab implementation of Diamant [START_REF] Charpentier | Diamanlab -An interactive Taylorbased continuation tool in MATLAB[END_REF] in comparison with these two classical first order continuation software. Then, the Elastica beam problem is used to exemplify our higher-order continuation framework on a classical partial differential equation (PDE) problem, relevant for numerical bifurcation analysis [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Doedel | Lecture notes on numerical analysis of nonlinear equations[END_REF][START_REF] Farrell | The computation of disconnected bifurcation diagrams[END_REF]] and structural mechanics.

The layout of the paper is as follows. Fundamentals of the higher-order continuation tool Diamant are briefly recalled in Section 2. The higher-order bifurcation analysis is presented in Section 3, then illustrated by several case studies in Section 4. Section 5 provides a summary.

Higher order continuation framework

Numerical continuation is a method devoted to the solution for the system of nonlinear parametric equations [START_REF] Allgower | Numerical continuation methods: an introduction[END_REF]. The solutions of such an under-determined problem are one-dimensional continua, referred to as solution branches, that may intersect at bifurcation points. The pseudo-arc length continuation algorithm [START_REF] Keller | Lectures on numerical methods in bifurcation problems[END_REF] is a classical answer to the calculation of solution branches. The branches U (a) = (u(a), λ(a)) are locally parameterized as

(U (a) -U (0)) * . ∂U ∂a -a = 0, ( 2
)
where * is the transpose operator and ∂U /∂a is the normalized direction vector. Under analyticity assumptions, the solution branches may be approximated as a collection of continuous truncated Taylor expansions [START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF][START_REF] Cochelin | Méthode asymptotique numérique[END_REF],

U (a) = K k=0 a k k! ∂ k u ∂a k (0), K k=0 a k k! ∂ k λ ∂a k (0) = K k=0 a k (u k , λ k ) = K k=0 a k U k , ∀a ∈ (0, a m ), (3)
where K is the truncation order of the series (typically K = 20), U k = (u k , λ k ) are the unknown Taylor coefficients of U (0) for k = 0, .., K and a m is the range of validity for the series. As in other higher-order solvers (see [START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF] and the references therein), this series is introduced into (1) to deduce a general iterative sequence of equations satisfied by the Taylor coefficients U k , that is

(R • U ) k = {R 1 }U k + {(R • U ) k |U k = 0} = 0, for k = 1, .., K. (4)
Therein, the Jacobian {R 1 } of R is constructed once. The higher-order term

{(R • U ) k |U k = 0}
is the Taylor coefficient at order k of (R • U ) evaluated with a null value for the unknown U k . All these derivatives may be computed using HOAD.

The solution for (1)-( 2) then relies on the sequence of linear systems

{R 1 } U * 1 .U k = -{(R • U ) k |U k = 0} δ 1k , for k = 1, .., K, (5)
where δ 1k is the Kronecker's delta. The left-hand side matrix of ( 5) is factorized once.

The range of validity a m is deduced from the remainder of the series as

a m = ε || {(R • U ) k |U k = 0} || 1/K , (6)
to provide large adaptive stepsize. The threshold ε is typically set to 10 -6 . This higherorder continuation method may often be used without Newton corrections.

Listing 1 Nonlinear solver for the computation of one branch of solutions % Given a solution point U 0, a tangent U t , a nonlinear function R Compute the Jacobian {R 1 } Form the left hand side matrix of (5) and factorize it Set

U 0 = U 0; U 1 = U t for k =2: K Evaluate R(U ) at order k using U k = 0 Solve (5) Set U k end 10
Compute a m following (6) % Room for HOAD simple bifurcation detection, see List. 4

The iterative sequence (5) is the basis for the general high level nonlinear solver known as Diamant. Its HOAD is similar to the HOAD described in other Taylor-based nonlinear solvers. Thanks to AD generality, this solver is able to deal with any user-defined analytical nonlinear residual problem satisfying [START_REF] Allgower | Numerical continuation methods: an introduction[END_REF]. The series information -Taylor coefficients, tangent to the branch and range of validity -are stored into checkpoints [START_REF] Charpentier | Checkpointing schemes or adjoint codes: Application to the meteorological model Meso-NH[END_REF][START_REF] Griewank | Algorithm 799: Revolve: An implementation of checkpoint for the reverse or adjoint mode of computational differentiation[END_REF] managed through the continuation driver described in List. 2. Singular point information discussed at Section 3 -fold and bifurcation loci, tangent to the second branch, series and related range of validity -is accounted in the checkpoint where relevant. This offers the opportunity to start a continuation at a bifurcation point in the direction of the second tangent, see List. 5. Previous Diamant versions did not provide a higher-order bifurcation analysis to improve the utility of the series near to a bifurcation point. This significant shortcoming is addressed by reworking the bifurcation strategies presented in [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Doedel | Numerical analysis and control of bifurcation problems (i) bifurcation in finite dimensions[END_REF] to propose a general higher-order AD implementation of bifurcation detection, location and branch switching.

Bifurcation analysis

Near to a singular point, a very small theoretical or numerical perturbation of the solution induces a qualitative change in the behavior of the residual function under study. Such a change can be monitored and classified from determinant calculations and/or derivative computations, see for instance [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory -Volume I[END_REF][START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]. Among the numerous kind of bifurcations, this paper consider fold points (or snap through in structural mechanics) and simple bifurcation points (or buckling) that are of special interest in PDE problems studied in structural engineering.

As an illustration, we consider the tutorial example of MatCont [START_REF] Kuznetsov | Tutorial II: One-parameter bifurcation analysis of equilibria with matcont[END_REF], where the variables u 1 and u 2 denote the components of the unknown state vector u, and α is the parameter to be varied. A trivial solution is O 1 = (u 1 , u 2 , α) = (0, 0, 0). The branch issued from O 1 satisfies α exp(u 1 ) -u 1 = 0 and u 1 = u 2 . It exhibits a simple fold at (1,1/ exp(1)) and a simple bifurcation point at (3,3,3/ exp(3)). The bifurcation diagram computed using the Matlab implementation of Diamant, namely Diamanlab, is plotted in Fig. 1. The branch issued from O 1 requires the computation of 5 checkpoints only to run from u 1 = 0 to u 1 = 8 with large continuation steps. The fold (0.999998, 0.999998, 0.367879), indicated with a green +, is identified from the series computed at the second checkpoint (a m = 2.1068 at checkpoint 2). The simple bifurcation point, indicated with a red disc, is located at (3.000001, 3.000001, 0.149361). As expected, the accuracy in the results agrees with the threshold of the method. The branching method (List. 5) allows to compute a series (with a range of validity a m = 0.813490), then the solution point O 2 on the second branch, far from the bifurcation point. The second branch is traveled forward and backward.

0 = -2u 1 + u 2 + α exp(u 1 ), 0 = u 1 -2u 2 + α exp(u 2 ). (7)
Fold and bifurcation detection as well as branch switching methods are described in subsections 3.1, 3.2 and 3.3, respectively.

Fold detection and location

As recalled in [START_REF] Doedel | Numerical analysis and control of bifurcation problems (i) bifurcation in finite dimensions[END_REF], a solution 1) is a simple fold if λ (a f ) = ∂λ ∂a (a f ) = 0, for some a f ∈ (0, a max ). In other words, the tangent to the solution branch at that point is vertical on a projected bifurcation diagram that plots λ along the horizontal axis. When the branch of solutions is approximated using series (3), the fold may be detected by monitoring the sign of λ (a) between successive checkpoints. This derivative is deduced from the series of λ in a straightforward manner as Note that the curvature λ (a f ) at the fold can be approximated from the series as

U f = U (a f ) = (u(a f ), λ(a f )) of (
λ (a) = K k=1 ka k-1 λ k (a). ( 8
) U 1 B 2 b U U 1 B 1 B 1 2 B (a) U 1 B 2 b U U 1 B 1 U(0,s) B 1 2 B * (b)
λ (a f ) = K k=2 k(k -1)(a f ) k-2 λ k (a f ), (9)
and allows to characterize U f as a simple quadratic fold when λ (a f ) = 0. The fold detection summarized in List. 3 takes place in the continuation driver (List. 2).

Listing 3 Fold detection and location carried out in the continuation driver

% Given a checkpoint point U computed using List. 1, that is % known series U (a) = K k=0 a k U k , known range a m and % known sign for λ (0) Evaluate U (a m ) = K k=1 a k m U k % Solution point 5 Evaluate λ (a m ) = K k=1 ka k-1 λ k (a m ) % Derivative of λ % Fold detection if ( sign (λ (0)) ∼= sign (λ (a m ))) % Fold location Find a f such that K k=1 k(a f ) k-1 λ k (a f ) = 0 %dichotomy 10 Evaluate λ (a f ) = K k=2 (k -1)k(a f ) k-1 λ k (a f )
Update the checkpoint with fold information end

Bifurcation detection and location

Bifurcation and branch switching is a more complex task. In the perfect case, Fig. 2.a, the two solution branches B 1 and B 2 intersect at the bifurcation point U b . Tangents to the branches at U b are denoted by U t1 1 and U t2 1 . Even a very small error, inherent to floatingpoint calculations, may induce the branch switching presented in Fig. 2.b. Bifurcation and stability analysis is thus of prime importance in computational and engineering sciences [START_REF] Allgower | Numerical continuation methods: an introduction[END_REF][START_REF] Doedel | Numerical analysis and control of bifurcation problems (i) bifurcation in finite dimensions[END_REF][START_REF] Keller | Lectures on numerical methods in bifurcation problems[END_REF]. Various methods exist to detect and to locate a bifurcation, to pass through it along the traveled solution branch or to switch from one branch to the other.

Unsurprisingly, such a singular behavior may be monitored by means of the derivatives of the nonlinear system, and notably using AD [START_REF] Dhooge | New features of the software matcont for bifurcation analysis of dynamical systems[END_REF][START_REF] Guckenheimer | Computing periodic orbits and their bifurcations with automatic differentiation[END_REF]. In this section, we take advantage of the HOAD computed Taylor series to detect a geometric series [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Van Dyke | Analysis and improvement of perturbation series[END_REF] that indicates a bifurcation, and to implement an efficient branch switching method at simple bifurcation points. The paper reworks using a HOAD formalism the method discussed in [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF] for quadratic formulations R(U ) = L 0 + L(U ) + Q(U, U ) (comprising a constant vector L 0 , a linear operator L and a quadratic operator Q) to extend it to the general nonlinear algebraic system (1).

Brief state of the art

Let U b be a simple bifurcation point on the solution branches B 1 and B 2 , Fig. 2.a. At point U b , the Jacobian {R b 1 } admits two null right vectors φ 1 and φ 2 and a left null vector ψ [START_REF] Doedel | Numerical analysis and control of bifurcation problems (i) bifurcation in finite dimensions[END_REF] that satisfy

{R b 1 }φ i = 0, for i = 1, 2, and ψ * .{R b 1 } = 0. ( 10
)
The null space and its orthogonal complement are denoted by

N ({R b 1 }) = span{φ 1 , φ 2 } and N ⊥ ({R b 1 }) = R n \ N ({R b 1 }), respectively. Let (α, β) ∈ R 2 and U b 1 be such that U b 1 = (αφ 1 + βφ 2 ) ∈ N ({R b 1 }). (11)
Differentiating (1) twice at point U b ,

{R b 1 }U b 2 + {R b 2 }U b 1 U b 1 = {R b 1 }U b 2 + {R b 2 }(αφ 1 + βφ 2 )(αφ 1 + βφ 2 ) = 0, (12)
and multiplying the result by ψ * yield

ψ * .{R b 2 }(α 2 φ 1 φ 1 + 2αβφ 1 φ 2 + β 2 φ 2 φ 2 ) = 0, (13)
and the so-called Algebraic Bifurcation Equation (ABE)

α 2 c 11 + 2αβc 12 + β 2 c 22 = 0, where c ij = ψ * .{R b 2 }φ i φ j for i = 1, 2. ( 14
)
The solutions of the ABE allow for the determination of the two tangents. The practical implementation reported in [START_REF] Doedel | Numerical analysis and control of bifurcation problems (i) bifurcation in finite dimensions[END_REF] monitors the determinant of {R 1 } along the traveled branch to detect a possible bifurcation between two successive continuation points. As demonstrated in [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF], a geometric series emerges between successive Taylor coefficients near to a bifurcation point. This particular event may be used for detection and branch switching. Four stages are necessary. Firstly, the Taylor coefficient sequence is monitored to detect a possible geometric series. Secondly, the geometric series is taken out to restore a clean series with an optimal range of validity along the traveled branch. Then, the bifurcation locus is determined. Thirdly, the tangent to the second branch is computed. Finally, an optimal (large) range of validity can be computed at the bifurcation point.

Evidence for a geometric series

Choosing φ 1 = U t1 1 and φ 2 = U t2 1 in (11) allows to write the first order representation of branches B 1 and B 2 as

U (α, β) = U b + αU t1 1 + βU t2 1 . (15)
In the perfect case, Fig. 2.a, the ABE simplifies into αβ = 0, and either α or β is null. This agrees with the particular choices done for φ 1 and φ 2 . The approximate solution U (α, β) theoretically lies on one or the other branch.

Numerical residual errors, Fig. 2.b, are usually not zero from a computer point of view. The ABE becomes αβ = µ where µ is a small real number. Using α = a -s and β = µ/(a -s) in ( 15) allows to write

U (a, s) = U b -sU t1 1 - µ s U t2 1 U (0, s) +aU t1 1 + µ s a a -s U t2 1 , (16)
where the shift s represents the distance from U (0, s) to U b . Equation ( 16) provides a first order representation of the perturbed branch B 1 evaluated at a point U 0 near to the simple bifurcation under study. Since µ/s is very small, the contribution of U t2

1 is negligible except when a is almost equal to s. In that case, the rational fraction a/(a -s) may be written as the geometric series k (a/s) k with common ratio a/s. This happens more easily for a small distance s and/or for a large residual error contribution µ at point U b . The interested reader is referred to [START_REF] Van Dyke | Analysis and improvement of perturbation series[END_REF] for the analysis of functional singularities through series.

The formation of a geometric series is monitored in the highest four Taylor coefficients of the last computed solution for [START_REF] Charpentier | Arbogast, higher-order AD for special functions with Modular C[END_REF]. The candidates as common ratio and scale factor are

σ = U * K U K U * K-1 U K , (17)
and the Taylor coefficient U K , respectively. A geometric series is detected when collinearity [START_REF] Van Dyke | Analysis and improvement of perturbation series[END_REF] and proportionality [START_REF] Farrell | The computation of disconnected bifurcation diagrams[END_REF] properties are satisfied

3 k=1 1 - U * K-k U K ||U K-k ||.||U K || 2 ≤ ε coll = 10 -4 , (18)
and

3 k=1 |(U K-k -σU K-k+1 ) * .U K-k | U * K-k .U K-k ≤ ε prop = 10 -3 . (19)
In practice, the geometric series is detected along the traveled branch either before or after the bifurcation point U b , depending on the distance s to U b and on the norm of the residual.

Cleaning of the current Taylor series and bifurcation location

The identified geometric series is taken off the Taylor coefficients to improve the utility of the series, that is to avoid (i) interactions between the two solution branches, (ii) a possible amplification of computational errors, (iii) small ranges of validity, or even (iv) an undesired branch switching as in Fig. 2.b.

and

N ⊥ ({R b 1 }), ψ * .({R b k+1 |U b k = 1, U b k+1 = 0}U b k + {R b k+1 |U b k = U b k+1 = 0}) = 0, (23)
to obtain the iterative sequence of linear systems satisfied by the unknowns U b k , k ≥ 2. The path equation contribution is

(U b 1 ) * .U b k = 0. ( 24
)
The solution for ( 22)-( 24) is managed assuming

U b k = α k U t1 1 + β k U t2 1 + Υ k , (25)
where the unknowns α k and β k are real parameters and Υ k ∈ N ⊥ ({R b 1 }). Firstly, the linear system for Υ k is deduced from the projection [START_REF] Govaerts | Matcont: Numerical bifurcation analysis toolbox in matlab[END_REF] and the path equation, that is

  {R b 1 } ψ (U t1 1 ) * 0 (U t2 1 ) * 0   .   Υ k γ   =   -{R b k |U b k = 0} 0 0   , (26)
where the (n + 2) × (n + 2) left-hand side matrix is nonsingular and γ is a Lagrange multiplier. Secondly, replacing U b 1 by U t1 1 and using ( 25) into the path equation ( 24) allow to deduce

(U t1 1 ) * .U b k = α k (U t1 1 ) * .U t1 1 + β k (U t1 1 ) * .U t2 1 = 0, (27)
that is

α k = -β k (U t1 1 ) * .U t2 1 (U t1 1 ) * .U t1 1 . (28)
Finally, equation [START_REF] Griewank | Algorithm 799: Revolve: An implementation of checkpoint for the reverse or adjoint mode of computational differentiation[END_REF] allows to determine β k and the final equation for U b k ,

U b k = β k U t2 1 - (U t1 1 ) * .U t2 1 (U t1 1 ) * .U t1 1 + Υ k . (29)
The optimal range of validity for the second series is deduced from the Taylor coefficient U b k following [START_REF] Charpentier | Sensitivity computations in higher order continuation methods[END_REF]. These different stages are implemented as in Listing 5.

Numerical results

This section first presents numerical results obtained for two tutorial examples from Mat-Cont and Auto. Numerical experiments are carried out using matcont5p3 [START_REF] Govaerts | Matcont: Numerical bifurcation analysis toolbox in matlab[END_REF], AUTO-07P [START_REF] Doedel | Collaborators[END_REF], and the Matlab implementation of Diamant called Diamanlab [START_REF] Charpentier | Diamanlab -An interactive Taylorbased continuation tool in MATLAB[END_REF]. These three software are usable for research purposes. Simulations are effected on commodity hardware comprising an Intel processor at a maximal frequency of 2.9 GHz.

October 21, 2019

Optimization Methods & Software HALversion˙OMS˙R2

Listing 5 Series computation at a bifurcation point, the transpose operator is here denoted by ' as in Matlab % Given the bifurcation point U = U b , the tangents U1t1 and U1t2, % and the nonlinear function R Compute the Jacobian {R 1 } and the left null vector ψ Factorize the left -hand side matrix of ( 26) Before proceeding, major practical differences in the implementations should be noted. Among them, Matcont and Auto solve ordinary differential equations while Diamanlab solves algebraic equations. Note that AUTO embeds a collocation method (spatial discretization). MatCont and Diamanlab deal with Matlab and propose interactive graphical interfaces, while Auto deals with Fortran codes and command lines. The use of software interfaces (ergonomy is out of scope here) and executable codes for matcont5p3 and AUTO-07P interfere with time measurements. Consequently, the software are mainly compared in terms of number of continuation steps.

Set U 0 = U (0); U 1 = U 1t1; % Order 2 Set V 2 = U 2 computed from (29) with β 2 = 0 and Υ 2 = 0 Project V 2 to compute pV 2 = ψ .{R 3 |U 2 = V 2 , U 3 = 0} Solve (26) to get Υ 2 Set U 2 = Υ 2 and evaluate W 3 = {R 3 |U 3 = 0} Project W 3 to get pW 3 Compute β 2 = -pW 3 /
The clamped-clamped Elastica beam [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Doedel | Lecture notes on numerical analysis of nonlinear equations[END_REF][START_REF] Farrell | The computation of disconnected bifurcation diagrams[END_REF] then illustrates the interest for higherorder continuation and bifurcation methods for PDE problems in structural mechanics.

Stationary solutions of Bratu-Gelfand problems

The Bratu-Gelfand problem,

∂u ∂t = ∆u + λ exp(u), t ≥ 0, x ∈ (0, 1), u(0) = u(1) = 0, ( 30
)
is a classical case study in continuation, see [START_REF] Doedel | Numerical analysis and control of bifurcation problems (i) bifurcation in finite dimensions[END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] for instance, and some numerical results have already been reported for the three software packages [START_REF] Charpentier | Sensitivity computations in higher order continuation methods[END_REF][START_REF] Eusebius | AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Eq February[END_REF][START_REF] Kuznetsov | Tutorial II: One-parameter bifurcation analysis of equilibria with matcont[END_REF].

Two internal point discretization

The tutorial [START_REF] Kuznetsov | Tutorial II: One-parameter bifurcation analysis of equilibria with matcont[END_REF] studies equilibrium solutions for a discretization of (30) by means of two internal points. In MatCont, the nonlinear system is implemented as

u 1 = -2u 1 + u 2 + α exp(u 1 ), u 2 = u 1 -2u 2 + α exp(u 2 ), ( 31
)
by using α = λ/9 and a symbolic differentiation for the first two derivatives. MatCont saves symbolically computed derivatives in a Matlab function. AD was not introduced in the ODE version of MatCont because it turned out to be slower than symbolic derivatives [START_REF] De | Numerical computation of normal form coefficients of ODEs in Matlab[END_REF]. Note that the (small) time necessary to precompute the model derivatives when using Matcont is not indicated and is not accounted at simulation time. The equilibrium solutions of (31) satisfy [START_REF] Charpentier | Diamanlab -An interactive Taylorbased continuation tool in MATLAB[END_REF] as well. Default initial (InitStepsize) and maximum stepsize (MaxStepsize) values are equal to 0.01 and 0.1. In comparison, assuming that the computed a m has a finite value, there is no limitation on the maximum range of validity for the series in Diamanlab. Solution branches plotted in Fig. 3.a show the fold (LP) and the branch point (BP). As MatCont solves the ODE system (31), it can also find out its Hopf bifurcation (H). The two methods, Figs. 1 and3.a, are equally accurate at computed points and singular points since they satisfy the same threshold.

For numerical experiments, MatCont is then run using various stepsizes to enable a comparison between the adaptive stepsize methods of MatCont and Diamanlab. Computed solution points are indicated with + in Figs. [START_REF] Charpentier | Checkpointing schemes or adjoint codes: Application to the meteorological model Meso-NH[END_REF].a and 3.b. Clearly, see Fig. 3.b, using a few larger stepsizes in a discrete method degrades the graphical representation. In comparison, a few series are sufficient to be obtain a very accurate plot, even near to singular points.

The numbers of continuation steps and checkpoints computed with MatCont and Diamant are reported in Tab. 1. Results indicate that, for an equivalent quality in the plots, the HOAD continuation can be carried out using less steps, less Jacobian calculations and less matrix factorizations, than a first order continuation method. This can be a valuable asset when the factorization of the left-hand side of ( 5) is more expensive than the computation of the higher order right-hand side terms. This is not the case in Matlab where performance and operator overloading (OO) look mutually exclusive. Time measurement is of difficult handling when dealing with interactive tools and so small

General case

On the one hand, the exp demo proposed in AUTO-07P [START_REF] Eusebius | AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Eq February[END_REF] computes the equilibrium solutions of (30) by means of the system of first order ODEs

   u = v, v = -λ exp(u), u(0) = u(1) = 0, (32)
where u and v are functions of the variable x only. These are discretized by means of a collocation method involving 10 mesh intervals and 4 Gauss points per interval. On the other hand, using Diamanlab, we consider the stationary discrete system comprising N equidistributed points

   R 1 (u, λ) = (-2u 1 + u 2 )/h 2 + λ exp(u 1 ), R i (u, λ) = (u i-1 -2u i + u i+1 )/h 2 + λ exp(u i ), for i = 2, .., N -1, R N (u, λ) = (u N -1 -2u N )/h 2 + λ exp(u N ), (33)
where h = 1/(N + 1) and u i , for i = 1, .., N , are the unknowns at discretization points.

The run of the exp demo is carried out from the trivial null solution as described in [START_REF] Eusebius | AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Eq February[END_REF]. Within Auto, the adaptive discretization method used 10 mesh intervals and 4 Gauss points per interval. Min and max stepsizes are set to 10 -3 and 2.10 -1 respectively. One hundred continuation steps allow to reach the last point, at label 3 in Fig 4(a). The fold at label 2 is reached after 33 continuation steps. Diamanlab is run using a finite difference discretization involving N = 39 equidistributed points to work with a similar number of unknowns. The fold is detected in the series of the second checkpoint. Bifurcation diagrams (λ, u max ) are plotted in Fig. 4 using Plaut and Matlab, respectively. Although computing six checkpoints is sufficient, see Tab. 2, any Matlab interpreted implementation can compete with a compiled code as far as computer time is considered.

Stationary solutions of an enzyme model

The two-compartment enzyme model (Kernevez,1980) implemented in the enz demo of AUTO-07P,

     s 1 = (s 0 -s 1 ) + (s 2 -s 1 ) - 100s 1 1 + s 1 + s 2 1 , s 2 = (s 0 -s 2 ) + (s 1 -s 2 ) - 100s 2 1 + s 2 + s 2 2 , (34)
is now chosen as a case study to compare the three software. Therein concentrations s 1 and s 2 are computed with respect to the variation of the concentration s 0 . Simulations are carried out using the three software to observe the number the continuation steps necessary to reach the singular points as labeled in Fig. 5(a) by using AUTO-07P. Surprisingly, the enz demo also makes use of an a priori non-necessary discretization (NTST=15 mesh intervals, NCOL= 4 Gauss collocation points). The maximum stepsize is set to 0.25 for both AUTO-07P and MatCont simulations.

Bifurcation diagram produced by Matcont and Diamanlab are plotted to observe the Hopf points detected by MatCont, Fig. 5(b) and the checkpoints computed by Diamanlab, Fig. 5(c). In this figure, the second branch is computed the point BP2 for a better rendering (a m = 2.81 between BP2 and O 2 ). Branching at BP1 yields a series with a smaller range of validity of a m = 0.092 because of neighboring singularities.

Conclusions raised from the first case study still apply, the series computed at the bifurcation point following List. 5 may have a large range of validity. Time measurements reveal that MatCont uses 0.4 s to compute the first 150 steps (just before the first fold), but it requires 1.5 s when it computes and plots the solution simultaneously. In comparison, Diamanlab computes and plots the first 8 checkpoints (outreaching the fold) in 1.57 s. The overhead due to OO is acceptable in an interactive continuation.

Elastica beam

We now consider the classical clamped-clamped inextensible elastic beam subject to a uniaxial compression [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Doedel | Lecture notes on numerical analysis of nonlinear equations[END_REF][START_REF] Farrell | The computation of disconnected bifurcation diagrams[END_REF], the bifurcation diagram of which contains many interconnected solution branches, Fig. 6. This case study is representative of large nonlinear models related to the static behavior of mechanical structures discretized by means of a finite element method. where (x(s), y(s)) is the position of the beam centerline, θ(s) is the rotation of the beam cross-section and m(s) is the bending momentum. Boundary conditions are x(0) = y(0) = θ(0) = 0 and m(0) = 2πC on the one side, and y(1) = θ(1) = 0 on the other side. The compressive load P varies while T and C are the lateral and torque reactions at the boundary. The ODE system is discretized on a regular mesh of N e elements with a cubic interpolation of the four unknown variables. An orthogonal collocation method at Gauss points is used to build an algebraic system that comprises about 3 × N e × 4 equations (55 in the present simulation). The truncation order, the threshold (6) and the Newton-Raphson threshold are set to K = 20, ε = 10 -6 and 2.10 -5 , respectively.

Simulation results

The simulation is carried out as follows. The initial point is located at (0.030,1.000) on the bifurcation diagram that plots the coordinate y(1/4) with respect to the variations of the load P . The initial tangent points downward. The lower part of the blue solution branch is obtained in 4 continuation steps. Checkpoints are indicated by crosses. The upper part is computed from the initial point by reversing the sign of the tangent at point (0.030, 1.000). Three more steps are computed. This blue curve corresponds to the first buckling branch and presents 3 bifurcation points indicated with circles and located at points (1.000,0.000) and (2.184,±0.195).

For visualization purposes, Diamanlab has been modified to change the line color to another one when a bifurcation point is selected from the graphical user interface as the beginning for a continuation along the second tangent. The red closed loop is plotted from the bifurcation located on the blue curve at point (2.184,0.195) in about 40 Taylor series computations. Bifurcations appear at (-1.400,±0.130). One observes that the distance between successive checkpoints varies and becomes smaller near to sharp folds. Bifurcations are detected and located on the fly without the need for small steps.

The green curve corresponds to the second buckling branch. It is plotted from the bifurcation located at (-1.400,0.130) on the red loop. About 30 checkpoints are computed and a new bifurcation is found. Former observations about steps size near to folds and bifurcations are confirmed. Finally the trivial compression branch is followed and plotted as a black line. The bifurcation diagram can be built from the trivial branch too. The complex bifurcation diagram is accurately and continuously handled with the proposed HOAD continuation and bifurcation framework. This approach is especially interesting for large structures involving numerous degrees of freedom since the continuation steps between successive checkpoints are large, notably near to the bifurcations. The number of Jacobian calculations and linear system solutions is thus limited.

Discussion

The use of higher order methods (continuation and bifurcation analysis) yields continuous solution branches that can be plotted as they are in a bifurcation diagram. It also reduces the number of checkpoints that has to be computed, including near to the singular points.

The efficiency in terms of computer time of this Matlab prototype of Diamant is obviously penalized by the use of operator overloading. In spite of this, the compromise between interactivity, generality and time remains actually acceptable for prototyping numerical methods and for educational purposes, notably in computational mechanics. In particular, we plane to interface Diamanlab with PDE codes involving plate or shell finite elements through the user-defined R.m method that implements the nonlinear function R of (1).

To that end, time performances can be optimized in several manners. First, sparse Jacobian calculations or, better, Jacobian calculation at assembly time can be implemented for PDE problems. Second, the HOAD computations in ( 5) and ( 26) could be improved to avoid re-computations by working on the HOAD library to preserve generality. Finally, HOAD and Diamant can be implemented in another language.

An attractive option will be to consider Arbogast [START_REF] Charpentier | Arbogast, higher-order AD for special functions with Modular C[END_REF], the toolbox for HOAD based on Modular C [START_REF] Clauss | Iterative Computations with Ordered Read-Write Locks[END_REF], as it provides an interface to differentiate C programs and to implement new differential operators from scratch. As Diamant solvers for continuation [START_REF] Charpentier | Arbogast, higher-order AD for special functions with Modular C[END_REF], bifurcation analysis [START_REF] Keller | Lectures on numerical methods in bifurcation problems[END_REF], homotopy [START_REF] Bilasse | A generic approach for the solution of nonlinear residual equations. Part II: Homotopy and complex nonlinear eigenvalue method[END_REF] and sensitivity analysis [START_REF] Charpentier | Sensitivity computations in higher order continuation methods[END_REF], are differential operators, they can be good candidates to benefit from Arbogast in terms of code and performance optimization.

Conclusions

Theoretical aspects of higher-order continuation and bifurcation analysis are presented in a HOAD formalism to emphasize the generality of this framework. The key point is in the series computations that enable for large adaptive continuation steps notably near to bifurcation points, and that limit the number of solutions/checkpoints to be computed while preserving the quality of bifurcation diagram plots. Strengths and current limitations are discussed on classical case studies from Auto and MatCont. Then, the Elastica beam example shows that complex bifurcation diagrams and simple bifurcations are handled in an accurate and robust manner. The proposed method is especially interesting for large structures involving numerous degrees of freedom because the number of Jacobian calculations and linear system solutions get smaller as continuation steps between successive checkpoints are large.

Future work is concerned with the solution of PDE problems through a finite element modeling (plate and shell elements) by using Diamanlab. The implementation in Arbogast will be then studied.

Figure 1 .

 1 Figure 1. Tutorial bifurcation diagram using Diamanalab.

Figure 2 .

 2 Figure 2. Bifurcation: (a) Perfect case, (b) Numerical case.

Figure 3 .

 3 Figure 3. Tutorial bifurcation diagram using MatCont. (a): Default stepsizes. (b): Large Stepsizes.

FoldFigure 4 .

 4 Figure 4. Bifurcation diagram for the equilibriums of (30). (a) AUTO-07P, (b) Diamanlab.

Figure 5 .

 5 Figure 5. Bifurcation diagram for the enz demo. (a) AUTO-07P, (b) MatCont, (c) Diamanlab from BP2.

Figure 6 .

 6 Figure 6. Bifurcation diagram for the clamped-clamped Elastica beam. The first two buckling branches (blue and green curves) are connected to the red closed loop at four bifurcation points.

Table 1 .

 1 Number of continuation steps/checkpoints computed with MatCont and Diamant for the two internal point discretization.

		MatCont's MaxStepsize Diamant
	InitStepsize 0.01 1	10	-
	MaxStepsize	0.1	1	10	-
	Branch 1	124	30 10	5
	Branch 2	210	52 13	17

Table 2 .

 2 Number of continuation steps/checkpoints and CPU times necessary to Auto and Diamanlab. systems. It results that Diamanlab's OO implementation needs 0.37 s to compute and to plot the first checkpoint, what is acceptable in an interactive continuation, where MatCont's needs 0.2 s to compute and plot the first 20 points.

		AUTO-07P		Diamanlab
		Label Nb steps time	Nb checkpoints time
	Fold	2	33		2
	Last point 3	100	0.044s 6	1.49 s

Table 3 .

 3 Enzyme model: Number of continuation steps/checkpoints necessary to Auto, MatCont and Diamanlab.

	Branch 1 Type Label AUTO-07P	matcont5p3 Diamanlab
	LP	2	141	152	8
	H	-	-	155	-
	BP1	3	145	157	10
	BP2	4	214	211	18
	H	-	-	224	-
	LP	5	244	243	19
	Branch 2 Type Label AUTO-07P	matcont5p3 Diamanlab
	H	-	-

The clean series U (a),

provides the perfect branch of solutions in the current tangent direction. Its range of validity is evaluated using ( 6) at order K -1. By construction, the bifurcation point satisfies U b = U (σ). An approximation of the tangent to the traveled branch at that point is

Determination of the second tangent

At the bifurcation point U b , the Taylor coefficient U t2 1 defining the second tangent belongs to the kernel N ({R b 1 }) and satisfies the ABE. Using U t2 1 = αU t1 1 + βφ 2 with φ 2 belonging to N ({R b 1 }) allows to compute U t2 1 from (13).

Branch switching

The series at the bifurcation point U b may be written as

where U b 1 is one of the two tangents and U b k , for k = 2, .., K, are unknown Taylor coefficients. This series is introduced into (1). Resulting equations are projected onto N ({R b 1 }),

Equations and discretization

Let s ∈ (0, 1) be the dimensionless curvilinear coordinate. The beam equilibriums satisfy

x (s) = cos(θ(s)), y (s) = sin(θ(s)), θ (s) = m(s), m (s) = -4π 2 (P sin(θ(s)) + T cos(θ(s)),