# Numerical study of the stability of solutions for the half-space Ginzburg-Landau model

Pierre del Castillo, Séraphin Mefire

► To cite this version:

Pierre del Castillo, Séraphin Mefire. Numerical study of the stability of solutions for the half-space Ginzburg-Landau model. Journal of Engineering Mathematics, 2011, 69 (4), pp.327-344. 10.1007/s10665-010-9380-1 . hal-01286184

HAL Id: hal-01286184

https://hal.science/hal-01286184

Submitted on 10 Mar 2016

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical Study of the Stability of Solutions for the Half-Space Ginzburg-Landau Model

P. Del Castillo* and S.M. Mefire†

**Abstract.** Based on a shooting alternative that allows us to numerically solve the one-dimensional system of Ginzburg-Landau in an unbounded domain, we perform here a numerical study of the stability of solutions of this system. This stability notion, from a physical point of view, means that each solution of the system is identified as stable when it minimizes the corresponding Ginzburg-Landau functional. As opposed to our previous paper, we are concerned here with a more general study since the weak and large regimes of the Ginzburg-Landau parameter are considered and the initial data are no longer subject to the de Gennes condition. We also numerically investigate certain conjectures regarding the superheating field.

**Keywords.** superconductivity, Ginzburg-Landau systems, initial value problems, shooting methods, Hermite finite elements, generalized eigenvalue problems

**MS Codes.** 34A34, 65L05, 65L06, 65N25, 65N30, 82D55

## 1 Introduction

The Ginzburg-Landau theory allows one to describe the states of a superconducting material in an exterior magnetic field $H_e$. When the material is a film delimited by two planes ($x = -d, x = d$), with the exterior field parallel to the surface, the determination of states from this theory involves a functional depending on the inner magnetic potential $A$ and on a real-valued wave function $f$ (cf. e.g. [1–5]). In fact, the states of the film are then characterized from global or local minima of this functional. When we are concerned with pairs $(f, A)$, with $f$ even and $A$ odd on $[-d, d]$, also called symmetric solutions, we thus consider the reduced Ginzburg-Landau functional,

$$ \varepsilon_d(f, A; h) = \int_0^d \left[ \frac{1}{2} f(x)^4 - f(x)^2 + \kappa^{-2} (f'(x))^2 + f(x)^2 A(x)^2 + (A'(x) - h)^2 \right] dx, $$

where $d$ is proportional to the thickness of the film and $h$ is proportional to $|H_e|$. In the definition of $\varepsilon_d$, it is assumed that the edge of the film is at $x = 0$. Also, $(f, A) \in (H^1((0, d)))^2$ with $A(d) = 0$. The so-called normal phase corresponds to the case where $f \equiv 0$ with $A' \equiv h$, whereas the superconducting phase corresponds to the context where $f \neq 0$. The Ginzburg-Landau parameter, $\kappa > 0$, characterizes the properties of the material. Its values determine the type of superconductor according to the type

<footnote>*LAMFA, CNRS, Université de Picardie, 33 r. Saint-Leu, 80039 Amiens Cedex 1, France. Email: Pierre.Delcastillo@u-picardie.fr</footnote>

<footnote>†IECN, CNRS, Université de Nancy-I, 54506 Vandœuvre-lès-Nancy Cedex, France. Email: Seraphin.Mefire@iecn.unancy.fr</footnote>
of phase transition which takes place between the normal phase and the superconducting phase: $\kappa$ small describes what is known as a type I superconductor and $\kappa$ large as a type II superconductor. The critical value of $\kappa$ that is usually given to separate type I and type II superconductors is $\kappa = \frac{1}{\sqrt{2}}$. For a type I superconductor, there is a critical value $h_c^*$ such that if $h < h_c^*$, the material is entirely superconducting and the magnetic field is expelled from the sample away from a boundary layer. This is called the Meissner effect. If $h > h_c^*$, superconductivity is destroyed and the material is in the normal state. Regarding applications, there have been proposals to use superheated type I superconductors as detectors for elementary particles, so that the sample acts as a superconducting “bubble chamber” (cf. e.g. [6]). The passage of a sufficiently energetic particle through the sample would initiate the transition to the normal state. For a type II superconductor, the phase transition is different and there are two critical values: $h_{c1}$ and $h_{c2}$. For $h < h_{c1}$, the exterior magnetic field is expelled from the sample and there is a Meissner effect as for type I superconductors. When $h$ increases above $h_{c1}$, superconductivity is not destroyed straight away, since the superconducting and the normal phases coexist in the form of filaments or vortices. As $h$ increases further, the vortices become more numerous until the critical value $h_{c2}$ is reached at which superconductivity is destroyed. For $h > h_{c2}$, the material is in the normal state. The way superconductivity is nucleated is highly dependent on $d$ and $\kappa$ (see e.g. [7, Chap. 3]). Further discussions of the development of the above model as well as for the non reduced model, and physical significances, are given in [2], [8–12].

When the width $d$ of the film is large (in the sense that $\kappa d$ is large), we are concerned with a functional, $\varepsilon_\infty$, formally introduced by setting $d = +\infty$ in the definition of $\varepsilon_d$ after a renormalization obtained by adding the term $(\frac{1}{2} - h^2)d$. In fact,

$$ \varepsilon_{\infty}(f, A; h) = \int_{0}^{+\infty} \left[ \kappa^{-2} f'(x)^2 + \frac{1}{2}(1 - f(x)^2)^2 + f(x)^2 A(x)^2 + A'(x)^2 \right] dx + 2hA(0), $$

where $(f, A) \in \mathcal{H}_\infty = \{(f, A); (1-f) \in H^1((0, +\infty)), A \in H^1((0, +\infty))\}$. The corresponding Ginzburg-Landau equations, expressing the necessary conditions for the existence of minima, are then

$$ (GL)_\infty \left\{ \begin{array}{ll} -\kappa^{-2}f'' - f + f^3 + f A^2 & = 0 & \text{on } (0, +\infty), \\ -A'' + Af^2 & = 0 & \text{on } (0, +\infty), \end{array} \right. \qquad (1) $$

with the initial boundary conditions

$$ f'(0) = 0, \quad A'(0) = h, \qquad (2) $$

and $(1 - f, A) \in (H^2((0, +\infty)))^2$. The problem (1) - (2) is called the half-space model and was studied in [3,4]. Following [13], any pair $(f, A)$ satisfying also (1) - (2) is such that $f > 0$ and is thus called a *positive solution*. On the other hand, $f \le 1$. With the choice of the variational space $\mathcal{H}_\infty$, the normal solutions are eliminated and a particular role is attributed to the positive solutions since it is then enforced that $f$ tends to 1, as $x$ tends to $+\infty$. For any fixed $\kappa$, we know (see [2]) that the set of $h$, such that there exists a solution of (1) - (2), is a bounded interval $[0, h^+]$. We then introduce the *superheating field*, $h^{sh}(\kappa)$, defined as the supremum of this interval. This critical field is very important
for many applications. For instance, measuring the superheating field provides one of the few methods for experimentally determining the characteristic of the superconductor (see [14, Chap. 16] for other properties and developments). The Ginzburg-Landau equations corresponding to $\epsilon_d$ constitute a system, denoted by $(GL)_d^s$, introduced similarly as $(GL)_\infty$ but on $(0, d)$ and with in addition to (2): $f'(d) = 0$, $A(d) = 0$. Let us recall that Bolley, Foucher and Helffer have shown in [1, 15] that the system (1) - (2) is close to $(GL)_d^s$ for $\kappa$ small, with $\kappa d$ large. More precisely, they have proved the existence of a local superheating field for $(GL)_d^s$ which is near to the global superheating field $h^{sh}(\kappa)$ for (1) - (2) when $\kappa d$ is large.

In [16], Bolley and Del Castillo have established the following existence result, regarding $(GL)_\infty$.

**Theorem 1.1** For all $\kappa > 0$ and for all $f_0 \in (0, 1]$, there exists a unique solution $(f, A)$ of (1) such that:

$$f(0) = f_0, \quad f'(0) = 0. \tag{3}$$

We can then introduce (see also [17]), for $\kappa > 0$,

$$\sigma_\kappa : f_0 \in (0, 1] \mapsto \sigma_\kappa(f_0) := A'(0),$$

called the *response map*. Hereafter, we denote a solution of (1) - (2) by $(f, A; h)$; the same notation will represent also a solution of (1) and (3), according to the response map.

A critical point $(f, A)$ of $\epsilon_\infty$ is said to be *stable* if the quadratic form $d^2\epsilon_\infty(f, A; h)$ is positive definite, and *unstable* if $d^2\epsilon_\infty(f, A; h)$ can take negative values.

In [18], we numerically studied the stability of solutions of (1) - (2), in the weak-$\kappa$ regime and in the context of a localization of these solutions. Precisely, we were concerned with pairs $(f, A)$ satisfying (1) - (2) and in accordance with the de Gennes condition (see [19]), $\kappa h^2 = \sqrt{2}f(0)^2(1 - f(0)^2)$, for $f(0) > \frac{1}{\sqrt{2}}$. The case $f(0) = \frac{1}{\sqrt{2}}$ was taken into account in the study of the stability when the intensity $h$, of the exterior magnetic field, was considered near the superheating field.

Based on a formal approach, namely the method of matched asymptotic expansions, Dolgert, Di Bartolo and Dorsey [20] introduce in the weak-$\kappa$ regime an expansion of $h^{sh}(\kappa)$ in powers of $\kappa^{1/2}$. It results from this paper that the numerical values of $h^{sh}(\kappa)$, obtained from a computation of solutions by a Newton-type method, agree quite well with the formal values given by the introduced asymptotic expansion. These authors have studied moreover the stability of solutions of (1) - (2) in the weak-$\kappa$ regime, with the same formal approach.

As opposed to [18], here we numerically study the stability of solutions of (1) - (2): in the weak and large-$\kappa$ regimes, in the context where the pairs $(f(0), h)$ are not subject to the de Gennes condition, and for $f(0) \in (0, 1)$.

In contrast with the approach of [20], we are here concerned with a numerical framework that leads us both to a determination of the superheating field and to a study of the stability of solutions in each regime of $\kappa$. This framework requires the combination of a shooting approach (introduced rigorously in [17]) with a semi-implicit method (which is A-stable) for the numerical computation of solutions, as well as Hermite element approximations in the stage of the stability study. The efficiency of this
framework is distinguished by the fact that small values of the initial datum $f_0 = f(0)$ can be here considered, contrary with an approach using the method of matched asymptotic expansions as in [20], where the formal pair, from which the expansion of $h^{sh}(\kappa)$ is determined, cannot be an approximate solution of (1) - (3) for such an initial value. The robustness of the present framework is noticed since it also allows us to study the stability of solutions between the two regimes of $\kappa$.

This paper is subdivided into five sections. In Section 2, we are concerned with some properties of the solutions of $(GL)_\infty$, and of $\sigma_\kappa$. After describing in Section 3 the numerical approximation of solutions of (1) - (2), based on a shooting alternative, we numerically inspect $\sigma_\kappa$ and other maps associated with $(GL)_\infty$, as well as certain conjectures regarding the superheating field. In Section 4, we numerically study the stability of the solutions of (1) - (2) in the weak and large-$\kappa$ regimes as well as between these regimes. Our numerical investigations lead us mainly to conjecture that the particular initial value $f_0$, from which $h^{sh}(\kappa)$ is determined, corresponds to the initial datum for which the change of stability of the associated solution occurs — independently of the values considered here for $\kappa$. Our conclusions are given in Section 5.

## 2 Preliminaries

Let us start by recalling some general properties (see [2, 21]) of the solutions of $(GL)_\infty$.

**Proposition 2.1** Let $(f, A)$ be a solution of (1) - (2). The pair $(f, A)$ is such that:

* $f$ is increasing on $[0, +\infty)$ and we have, for all $x \in [0, +\infty)$,

$$0 < f(x) \le \tanh(\frac{\kappa x}{\sqrt{2}} + x_0), \quad \tanh(x_0) = f(0), \qquad (4)$$

* $A$ is strictly negative and increasing on $[0, +\infty)$; we have

$$0 \le A'(x) \le h, \quad \forall x \in [0, +\infty),$$

* the following energy conservation law,

$$\kappa^{-2} f'(x)^2 + A'(x)^2 = A(x)^2 f(x)^2 + \frac{(1 - f(x)^2)^2}{2}, \quad \forall x \in [0, +\infty), \qquad (5)$$

is satisfied.

The following statement recalls the properties of the response map (see [13, 17], [22–24]).

**Proposition 2.2** We have:

$$\sigma_{\kappa} \in C^0((0, 1]), \qquad (6)$$

$$\lim_{f_0 \to 0} \sigma_{\kappa}(f_0) = \frac{1}{\sqrt{2}}, \quad \sigma_{\kappa}(1) = 0.$$
In the weak-κ regime, we have

$$ \max_{f_0 \in (0,1]} \sigma_{\kappa}(f_0) = 2^{-\frac{3}{4}} \kappa^{-\frac{1}{2}} \left( 1 + \frac{15}{32} \sqrt{2\kappa} + o(\kappa) \right). \quad (7) $$

In the large-κ regime, we have

$$ \lim_{\kappa \to +\infty} \max_{f_0 \in (0,1]} \sigma_{\kappa}(f_0) = \frac{1}{\sqrt{2}}. $$

**Remark 2.1** Considering (5) at the point $x = 0$, and according to the definition of $\sigma_{\kappa}$, it results that

$$ \sigma_{\kappa}(f_0)^2 = A_0^2 f_0^2 + \frac{(1 - f_0^2)^2}{2}, \quad \text{where } A_0 := A(0). \quad (8) $$

From Theorem 1.1, we can introduce the following maps

$$ \mathcal{F}: (x, f_0) \in \mathbb{R}^+ \times (0,1] \longmapsto f(x; f_0), \quad \mathcal{A}: (x, f_0) \in \mathbb{R}^+ \times (0,1] \longmapsto A(x; f_0), \quad (9) $$

where $f(x; f_0)$ and $A(x; f_0)$ represent the values of the solution $(f, A; h)$ at the point $x$, and for $f_0 \in (0,1]$ given. Based on (6), we can show that these maps have the following properties.

**Proposition 2.3** We have

$$ \mathcal{F} \in C^0(\mathbb{R}^+ \times (0,1]) \quad \text{and} \quad \mathcal{A} \in C^0(\mathbb{R}^+ \times (0,1]). $$

*Proof.* Let us introduce $G \in C^\infty(\mathbb{R}^4, \mathbb{R}^4)$, defined for $v = (v_1, v_2, v_3, v_4)^T \in \mathbb{R}^4$, as follows:

$$ G(v) = (v_2, \kappa^2(-v_1 + v_1^3 + v_3^2 v_1), v_4, v_1^2 v_3)^T, $$

where the superscript “T” denotes the transpose. For $f_0 \in (0,1]$, let us consider $(f, A; h)$ the solution of (1) and (3). By representing respectively by $f'(x; f_0)$, $A'(x; f_0)$ the values of $f'$, $A'$, at the point $x$, for $f_0$ then fixed, we notice that $Y(x; f_0) := (f(x; f_0), f'(x; f_0), A(x; f_0), A'(x; f_0))^T$, with $x \in \mathbb{R}^+$, is the solution of:

$$ \left\{ \begin{array}{l} \displaystyle \frac{dY}{dx} = G(Y) \quad \text{on } [0, \infty), \\ Y(0; f_0) = (f_0, 0, A_0(f_0), \sigma_{\kappa}(f_0))^T, \end{array} \right. \quad (10) $$

where, from Proposition 2.1 and (8), $A_0(f_0) := -\frac{\sqrt{\sigma_{\kappa}(f_0)^2 - \frac{(1-f_0^2)^2}{2}}}{f_0}$. Let us consider the map $\phi$ defined on $(0,1]$ as follows:

$$ \phi(f_0) := (f_0, 0, A_0(f_0), \sigma_{\kappa}(f_0))^T. \quad (11) $$

From (6) and (11), it results that $\phi$ is continuous on $(0,1]$. Let us fix $(x_1, f_1) \in \mathbb{R}^+ \times (0,1)$, and set $w := (f_1, 0, A_0(f_1), \sigma_{\kappa}(f_1))^T$, $K := [0, \alpha]$, where $\alpha > x_1$. Also, let us set, for $\eta > 0$, $B(w, \eta) := \{x \in \mathbb{R}^4; \|x - w\|_2 \le \eta\}$, where $\|x\|_2 := (\sum_{i=1}^4 x_i^2)^{\frac{1}{2}}$ with $x = (x_1, x_2, x_3, x_4)^T$. As $G \in C^\infty(\mathbb{R}^4, \mathbb{R}^4)$, from standard arguments, we can define, for some $\eta_\alpha > 0$, $\tilde{Y}: (x, u) \in K \times B(w, \eta_\alpha) \mapsto \tilde{Y}(x; u)$ such that:

$$ \left\{ \begin{array}{l} \displaystyle \frac{d\tilde{Y}}{dx} = G(\tilde{Y}) \quad \text{on } [0, \alpha], \\ \tilde{Y}(0, u) = u. \end{array} \right. \quad (12) $$
Moreover, this map is continuous on $K \times B(w, \eta_\alpha)$. On the other hand, from the continuity of $\phi$ on $(0, 1]$, there exists $\rho_\alpha > 0$ such that, for all $f_0 \in (f_1 - \rho_\alpha, f_1 + \rho_\alpha)$, we have $\phi(f_0) \in B(w, \eta_\alpha)$. We then deduce that the map $(x, f_0) \mapsto \tilde{Y}(x, \phi(f_0))$, defined on $K \times (f_1 - \rho_\alpha, f_1 + \rho_\alpha)$, is continuous at $(x_1, f_1)$. From (10) - (12), $Y(x; f_0) = \tilde{Y}(x, \phi(f_0))$ for $(x, f_0) \in K \times (f_1 - \rho_\alpha, f_1 + \rho_\alpha)$. Hence the maps defined in (9) are continuous at $(x_1, f_1)$. Similarly, we can obtain the continuity of these maps at $(x_1, 1)$. $\square$

# 3 Numerical inspections of certain maps and conjectures

In this section, we numerically study the behavior of the maps $\mathcal{F}$, $\mathcal{A}$ introduced in (9), and inspect certain conjectures regarding the superheating field. We are thus concerned at first with the numerical approximation of the solutions of $(GL)_\infty$ associated with the initial conditions in (2).

## 3.1 Shooting alternative and numerical approximations

A way to overcome the numerical difficulties due to the conditions at infinity consists of approximating the solutions of (1) - (2) from computations involving an associated initial value problem. Such a problem is obtained by replacing the conditions at $+\infty$ involving $f$ and $A$, by suitable initial conditions for $f$ and $A$. It consists of finding $f$ and $A$, for $(f_0, A_0; h)$ given, such that:

$$ (GL)_{ivp} \quad \left\{ \begin{array}{lclcl} -\kappa^{-2} f'' - f + f^3 + f A^2 & = & 0 & \text{on } (0, T_{max}), \\ -A'' + A f^2 & = & 0 & \text{on } (0, T_{max}), \\ f'(0) & = & 0, & A'(0) & = & h, \\ f(0) & = & f_0, & A(0) & = & A_0, \end{array} \right. $$

where $I_{max} := [0, T_{max})$ is the maximal interval of existence. As any pair $(f, A)$, solution of (1) - (2), satisfies the energy conservation law (5), we assume that

$$ h^2 = \frac{1}{2}(1 - f_0^2)^2 + A_0^2 f_0^2, \quad A_0 < 0, \quad f_0 \in (0, 1). \quad (13) $$

As indicated also in Proposition 2.1, the solutions of (1) - (2) are defined on $(0, +\infty)$ and bounded. This is not the case for most of the solutions of $(GL)_{ivp}$. In order to select the trajectories of $(GL)_{ivp}$ that correspond to those of (1) - (2) and (13), we make use of a shooting alternative, suitable in all the regimes of $\kappa$. This alternative is introduced from the definition of two disjoint open sets $V_1(f_0) := (\sigma_\kappa(f_0), +\infty)$ and $V_2(f_0) := (b_0(f_0), \sigma_\kappa(f_0))$, where $b_0(f_0) := 2^{-\frac{1}{2}}(1 - f_0^4)^{\frac{1}{2}}$, and is formulated as follows (see Bolley and Helffer [17]):

**Theorem 3.1 i).** If $h \in V_1(f_0)$, the solution $(f, A)$ satisfying $(GL)_{ivp}$ and (13) is such that $f$ is increasing and becomes larger than 1. Moreover, the function $A$ is strictly negative on $I_{max}$ and becomes decreasing.

ii). If $h \in V_2(f_0)$, the solution $(f, A)$ satisfying $(GL)_{ivp}$ and (13) has the property that there exist $0 < x_0 < z_0 < T_{max}$ such that $f < 1$ on $[0, z_0]$ and $f' < 0$ on $(x_0, z_0]$. Moreover, $A$ is increasing on $I_{max}$ and becomes greater than 0.
Based on this alternative, we numerically determine the particular value $\sigma_{\kappa}(f_0)$ of $h$ which leads to a solution of (1) - (2), for $f_0$ and $A_0$ considered as in (13). More precisely, after fixing $f_0 \in (0, 1)$ at first, for each $h$, we take $A_0$ in accordance with (13). We start at $x = 0$ with the initial data of $(GL)_{ivp}$, and compute the functions $f, f', A$ and $A'$ by iteration on an increasing sequence $(x_k)_k$, $x_k = \sum_{i=0}^k \zeta_i$, from a Runge-Kutta like method (with a variable step of size $\zeta_i$). By testing on such a sequence the criteria on $f, f', A$ and $A'$ provided from the contexts $i$) and $ii$) of Theorem 3.1, we can deduce whether the considered $h$ is above or below the value $\sigma_{\kappa}(f_0)$. We mention here that our simulations reveal that only criteria on $f, f'$ can be involved in the large-$\kappa$ regime and only criteria on $A, A'$ can be enforced in the weak-$\kappa$ regime. Concretely, for $h_1$ (large) in $V_1(f_0)$, and $h_2$ (small) in $V_2(f_0)$, we define $h_3 := \frac{1}{2}(h_1 + h_2)$. Following the context where $h_3 \in V_1(f_0)$, or $h_3 \in V_2(f_0)$, we define $h_4 := \frac{1}{2}(h_2 + h_3)$ or $h_4 := \frac{1}{2}(h_1 + h_3)$. This dichotomy process leads to two values $h_n$ and $h_{n+1}$, $h_n < h_{n+1}$, for some $n$, such that: $\sigma_{\kappa}(f_0) \in [h_n, h_{n+1}]$.

Before presenting the numerical method used for solving $(GL)_{ivp}$, let us reformulate this system as follows:

$$ \left\{ \begin{array}{ll} Y'(x) = F(Y(x), x) & \text{for } x \in (0, T_{\text{max}}), \\ Y(0) = Y_0, & \end{array} \right. \qquad (14) $$

where $Y(x) = (f(x), f'(x), A(x), A'(x))^T$, $Y_0 = (f_0, 0, A_0, h)^T \in \mathbb{R}^4$, and $F$ is a $\mathbb{R}^4$-valued function, defined on $\mathbb{R}^4 \times (0, T_{\text{max}})$. For $x \in (0, T_{\text{max}})$,

$$ F(Y(x), x) = (f'(x), \kappa^2(-f(x) + f(x)^3 + A(x)^2f(x)), A'(x), f(x)^2A(x))^T. $$

Each solution $(f, A)$ of (1) - (2) will be determined from a solution $Y$ of (14), when $Y = (Y_1, Y_2, Y_3, Y_4)^T$ is defined on $[0, +\infty)$ for $f_0, A_0, h$ considered as in (13), with $\lim_{x \to +\infty} Y_1(x) = 1$ and $\lim_{x \to +\infty} Y_3(x) = 0$. As in [25], in the context of a numerical study of the superheating field in the weak-$\kappa$ regime, we use the semi-implicit Runge-Kutta method of order 3, with two intermediate steps, described below. This is an A-stable method (see [26, Chap. 5]), which appears efficient for solving stiff problems. We are here concerned of course with an initial value problem which is very sensitive to small variations of the initial data, since we are looking for trajectories tending to the unstable stationary solution $(1, 0, 0, 0)^T$ as $x$ tends to $\infty$.

The numerical approximation of the solution $Y$ is achieved on an ordered sequence of points $x_i$, $i = 1, \dots, n$, of $(0, T_{\text{max}})$. Denoting by $Y_i$ an approximation of $Y(x_i) = (f(x_i), f'(x_i), A(x_i), A'(x_i))^T$ at $x_i$, $1 \le i \le n$ fixed, and by $H_i$ the variable steplength of the mentioned sequence, we determine the approximation $Y_{i+1}$ of $Y(x_{i+1})$ by solving the following system:

$$ \left\{ \begin{array}{l} Y_{i,1} = Y_i + \tau H_i F(Y_{i,1}, x_{i,1}), \\ Y_{i,2} = Y_i + \tau H_i F(Y_{i,2}, x_{i,2}) - \frac{1}{\sqrt{3}} H_i F(Y_{i,1}, x_{i,1}), \\ Y_{i+1} = Y_i + \frac{1}{2} H_i [F(Y_{i,1}, x_{i,1}) + F(Y_{i,2}, x_{i,2})], \end{array} \right. \qquad (15) $$
where

$$ \left\{ \begin{array}{l} \tau = \frac{1}{2}(1 + \frac{1}{\sqrt{3}}), \\ x_{i,1} = x_i + \tau H_i, \\ x_{i,2} = x_i + (1 - \tau) H_i. \end{array} \right. \qquad (16) $$

The parameter $H_i$ is determined at each step from comparative calculations involving two different steplengths and a local a posteriori estimate (see e.g. [27, Chap. 11] for such calculations). On the other hand, we use a Newton-type method for solving (15).

The present approach for approximating a solution of (1) - (2) is then based on the combination of the shooting alternative with this numerical method.

## 3.2 Numerical inspections

From the computation of each pair $(f, A)$ described as above, we numerically inspect here at first the behavior of $\mathcal{F}(x, f_0)$ and of $A(x, f_0)$, with respect to $f_0 \in (0, 1)$ as well as to $x \in \mathbb{R}^+$, in the weak and large-$\kappa$ regimes.

Figures 1 - 2 concern the weak-$\kappa$ regime. For each of the fixed values of $f_0$, we observe that $\mathcal{F}$ and $A$ increase with respect to $x$. This monotonicity of $\mathcal{F}$ and of $A$ were expected from Proposition 2.1. It occurs here that the behavior of $A$ seems to be more influenced by the values of $\kappa$, contrary to the behavior of $\mathcal{F}$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Behavior of $\mathcal{F}$ with respect to $x$, and for $f_0 = \frac{i}{20}\frac{1}{\sqrt{2}}$; $i = 1, 5, 10, 20$. At left $\kappa = 0.1$, and in the middle $\kappa = 0.05$. At right, representation (—) of $f_{lim}(x) = \tanh(\frac{\kappa x}{\sqrt{2}} + x_0)$, introduced in (4), and behavior (\dots) of $\mathcal{F}$ with respect to $x$, for $f_0 = \frac{i}{20}\frac{1}{\sqrt{2}}$, $i = 5$, $\kappa = 0.05$.</figcaption></figure>

In Figures 3 - 4, we represent the results obtained in the large-$\kappa$ regime. The same type of monotonicity with respect to $x$ and regarding $\mathcal{F}$, or $A$, is observed as in the weak-$\kappa$ regime.

Results similar to those represented in Figures 1 - 4 have been also obtained from simulations with $\kappa = 0.01, 0.001, 10, 25$. These results suggest moreover that, independently of $x$ and of the regime of
<figure><img src="image_3.png" /><figcaption>Figure 2: Behavior of $\mathcal{A}$ with respect to $x$, and for $f_0 = \frac{i}{20\sqrt{2}}; i = 1, 5, 10, 20$. At left $\kappa = 0.1$, and at right $\kappa = 0.05$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Behavior of $\mathcal{F}$ at left, and of $\mathcal{A}$ at right, with respect to $x$, for $f_0 = \frac{i}{20\sqrt{2}}; i = 1, 5, 10, 20$, and $\kappa = 1.3$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Behavior of $\mathcal{F}$ at left, and of $\mathcal{A}$ at right, with respect to $x$, for $f_0 = \frac{1}{\sqrt{2}} + \frac{i}{20}(1 - \frac{1}{\sqrt{2}})$; $i = 3, 7, 11, 15$, and $\kappa = 1.3$. In the middle, representation (—) of $f_{lim}(x) = \tanh(\frac{\kappa x}{\sqrt{2}} + x_0)$, introduced in (4), and behavior (\dots) of $\mathcal{F}$ with respect to $x$, for $f_0 = \frac{1}{\sqrt{2}} + \frac{i}{20}(1 - \frac{1}{\sqrt{2}})$, $i = 7$, $\kappa = 1.3$.</figcaption></figure>
$\kappa$, the maps $\mathcal{F}$ and $\mathcal{A}$ are increasing with respect to the considered values of $f_0$.

Let us now inspect the behavior of the response map $\sigma_\kappa$. The results in the weak-$\kappa$ regime are represented in Figure 5 whereas those obtained in the large-$\kappa$ regime are shown by Figure 6.

In the simulations, the sequence of values considered for $f_0$ is enriched, around its element providing a coarse evaluation of the superheating field, in order to determine with more accuracy the particular initial value that is associated with the superheating field.

<figure><img src="image_6.png" /><figcaption>Figure 5: Behavior of $\sigma_\kappa$ with respect to values of $f_0$ considered in $(10^{-4}, 1)$; $\kappa = 0.1, 0.3, 0.4$ at left, $\kappa = 0.01, 0.03, 0.04, 0.05$ in the middle, and $\kappa = 0.001, 0.003, 0.004$ at right.</figcaption></figure>

In each of the regimes of $\kappa$, it appears from Figures 5-6 that, on $(0, 1)$, $\sigma_\kappa$ admits a unique maximum: the superheating field $h^{sh}(\kappa)$. Moreover, the obtained results suggest that (1) - (2) has a unique solution if $h \le \frac{1}{\sqrt{2}}$, and two solutions if $\frac{1}{\sqrt{2}} < h < h^{sh}(\kappa)$.

<figure><img src="image_7.png" /><figcaption>Figure 6: Behavior of $\sigma_\kappa$ with respect to values of $f_0$ considered in $[10^{-3}, 1)$; $\kappa = 1.3, 3, 5$ at left, $\kappa = 9, 11, 15$ in the middle. Also, for $\kappa = 9, 11, 15$, a zoom detailing the behavior of $\sigma_\kappa$, with respect to small values of $f_0$, is performed at right.</figcaption></figure>

Before inspecting the conjectures regarding $h^{sh}$, let us examine how our numerical computations fit with the evaluation of $h^{sh}$ based on the Parr formula [28], which is recalled in (7) and rigorously proved in [23, 24]. In this way, we reproduce in the second column of Table 1 the values provided by the first two terms of the asymptotic expansion in (7), and represented here by $h_{Parr}^{sh} := 2^{-\frac{3}{4}}\kappa^{-\frac{1}{2}}(1 + \frac{15}{32}\sqrt{2}\kappa)$.
<table><caption>Table 1: Values of the superheating field, in the weak-κ regime, resulting from the Parr formula and numerical computations.</caption><thead><tr><th>κ</th><th>h<sup>sh</sup><sub>Parr</sub></th><th>h<sup>sh</sup><sub>num</sub></th><th>C<sup>num</sup><sub>Parr</sub></th></tr></thead><tbody><tr><td>0.4</td><td>1.189445</td><td>1.158998</td><td>0.581947</td></tr><tr><td>0.3</td><td>1.301488</td><td>1.280062</td><td>0.597122</td></tr><tr><td>0.1</td><td>2.004949</td><td>1.999854</td><td>0.635815</td></tr><tr><td>0.05</td><td>2.747287</td><td>2.745351</td><td>0.648350</td></tr><tr><td>0.04</td><td>3.051851</td><td>3.050444</td><td>0.651074</td></tr><tr><td>0.03</td><td>3.501217</td><td>3.500288</td><td>0.653887</td></tr><tr><td>0.01</td><td>5.985452</td><td>5.985267</td><td>0.659791</td></tr><tr><td>0.004</td><td>9.426437</td><td>9.426389</td><td>0.661629</td></tr><tr><td>0.003</td><td>10.877515</td><td>10.877484</td><td>0.661941</td></tr></tbody></table>

The third column of this table contains the numerical values of the superheating field, resulting from our computations, and represented by $h_{num}^{sh}$. Finally, the fourth column of the table contains the values of $C_{Parr}^{num} := \left( \frac{h_{num}^{sh}}{2^{-\frac{3}{4}}\kappa^{-\frac{1}{2}}} - 1 \right) \frac{1}{\kappa}$. In fact, when $\kappa$ becomes small enough, the comparison of $C_{Parr}^{num}$ with the constant $\frac{15}{32}\sqrt{2} \approx 0.662912$ appearing in (7) indicates that results fit well with the Parr formula, and validates in a certain sense our computations. For a more detailed numerical inspection of the Parr formula, we refer to [25].

The particular value of $f_0$, distinguished here by $f_0^{sh} := \frac{1}{\sqrt{2}} - \frac{7}{32}\kappa$, then associated (through $\sigma_\kappa$) with $h_{Parr}^{sh}$, derives from formal computations when $\kappa$ is small enough (see [20, 29]). We reproduce, with respect to $\kappa$, the values $f_0^{sh}(\kappa)$ in the second column of Table 2. The third column of this table contains the values of $f_0$ deriving from numerical computations, distinguished by $f_{0,num}^{sh}$, and that were associated with $h_{num}^{sh}$. It appears (see e.g. Figure 5) that $\sigma_\kappa$ is increasing on $(0, f_{0,num}^{sh}]$ and decreasing on $[f_{0,num}^{sh}, 1)$. The results of Tables 1 - 2 suggest in the weak-κ regime that there exist $\kappa_0$ and a map $\theta$, defined on $[0, \kappa_0]$, satisfying $\theta(\kappa) = \frac{1}{\sqrt{2}} - \frac{7}{32}\kappa + o(\kappa)$, such that for any $\kappa \le \kappa_0$, we have: $\sigma_\kappa(\theta(\kappa)) = \max_{f_0 \in [0,1]} \sigma_\kappa(f_0)$.

We are now interested in the numerical inspection of two conjectures: an asymptotic expansion of $h^{sh}$ proposed by Chapman [30], and an estimate (cf. [22]) that involves also $h^{sh}$.

The expansion of the superheating field conjectured in [30] concerns the large-κ regime and is expressed as follows:

$$h^{sh}(\kappa) = \max_{f_0 \in (0,1]} \sigma_{\kappa}(f_0) = \frac{1}{\sqrt{2}} + C\kappa^{-\frac{4}{3}} + o(\kappa^{-\frac{4}{3}}), \quad (17)$$

where $C \approx 0.33$.

Let us mention that the existence of a family of solutions of the second Painlevé equation which arose in the approach of [30] is investigated in [31], as well as some properties associated with the
<table><caption>Table 2: Values of $f_0^{sh}$ and of $f_{0,num}^{sh}$ in the weak-$\kappa$ regime ($f_{0,num}^{sh}$ being associated with $h_{num}^{sh}$ presented in Table 1).</caption><thead><tr><th>κ</th><th>f<sub>0</sub><sup>sh</sup></th><th>f<sub>0,num</sub><sup>sh</sup></th></tr></thead><tbody><tr><td>0.4</td><td>0.619606</td><td>0.642949</td></tr><tr><td>0.3</td><td>0.641481</td><td>0.656085</td></tr><tr><td>0.1</td><td>0.685231</td><td>0.687395</td></tr><tr><td>0.05</td><td>0.696169</td><td>0.696770</td></tr><tr><td>0.04</td><td>0.698356</td><td>0.698749</td></tr><tr><td>0.03</td><td>0.700544</td><td>0.700771</td></tr><tr><td>0.01</td><td>0.704919</td><td>0.704947</td></tr><tr><td>0.004</td><td>0.706231</td><td>0.706236</td></tr><tr><td>0.003</td><td>0.706450</td><td>0.706453</td></tr></tbody></table>

existence of the superheating field.

In order to numerically inspect formula (17), we reproduce in the second column of Table 3 the values of $h_{Chap}^{sh} := \frac{1}{\sqrt{2}} + 0.33\kappa^{-\frac{4}{3}}$ as well as in the fourth column, the values of $C_{Chap}^{num} := (h_{num}^{sh} - \frac{1}{\sqrt{2}})\kappa^{\frac{4}{3}}$, where the numerical values of the superheating field resulting from our computations, and also represented here by $h_{num}^{sh}$, are contained in the third column of this table.

The comparison of values of $C_{Chap}^{num}$ with the constant 0.33 suggests that, for $\kappa$ large enough, the formula proposed by Chapman appears to be here numerically satisfied.

The numerical results presented in Tables 1 and 3, as well as those obtained with the intermediate values, $\kappa = 0.6, 0.7, 0.707, 0.8, 0.9$, that are respectively $h_{num}^{sh} = 1.023937, 0.982080, 0.979540, 0.949568, 0.923575$, suggest that the following estimate involving the superheating field, and conjectured in [22],

$$h^{sh}(\kappa) > \frac{1}{\sqrt{2}}, \quad \forall \kappa > 0$$

appears also to be numerically satisfied.

By distinguishing also by $f_{0,num}^{sh}$, in Table 4 and for the large-$\kappa$ regime, the particular values of $f_0$ deriving from numerical computations and associated with $h_{num}^{sh}$ presented in Table 3, it appears that these values become (very slowly) smaller when $\kappa$ increases. Also, it results (see e.g. Figure 6) that $\sigma_\kappa$ is increasing on $(0, f_{0,num}^{sh}]$ and decreasing on $[f_{0,num}^{sh}, 1)$.

The results of Tables 3 - 4 and of Figure 6 lead us to conjecture in the large-$\kappa$ regime that there exists a map $\tilde{\theta}$, defined on $(0, +\infty)$, with $\tilde{\theta}(\kappa)$ decreasing when $\kappa$ becomes large, such that $\sigma_\kappa(\tilde{\theta}(\kappa)) = \max_{f_0 \in [0,1]} \sigma_\kappa(f_0)$.
<table><caption>Table 3: Values of the superheating field, in the large-$\kappa$ regime, resulting from numerical computations.</caption><thead><tr><th>κ</th><th>h<sup>sh</sup><sub>Chap</sub></th><th>h<sup>sh</sup><sub>num</sub></th><th>C<sup>num</sup><sub>Chap</sub></th></tr></thead><tbody><tr><td>1.3</td><td>0.939695</td><td>0.856835</td><td>0.212436</td></tr><tr><td>3</td><td>0.783376</td><td>0.766994</td><td>0.259116</td></tr><tr><td>5</td><td>0.745703</td><td>0.739820</td><td>0.279694</td></tr><tr><td>9</td><td>0.724734</td><td>0.722936</td><td>0.296334</td></tr><tr><td>10</td><td>0.722424</td><td>0.720968</td><td>0.298630</td></tr><tr><td>11</td><td>0.720596</td><td>0.719393</td><td>0.300567</td></tr><tr><td>15</td><td>0.716027</td><td>0.715378</td><td>0.305978</td></tr><tr><td>100</td><td>0.707817</td><td>0.707799</td><td>0.321299</td></tr></tbody></table>

<table><caption>Table 4: Values of $f_{0,num}^{sh}$ in the large-$\kappa$ regime (associated with $h_{num}^{sh}$ presented in Table 3).</caption><thead><tr><th>κ</th><th>1.3</th><th>3</th><th>5</th><th>9</th><th>10</th><th>11</th><th>15</th></tr></thead><tbody><tr><th>f<sup>sh</sup><sub>0,num</sub></th><td>0.561301</td><td>0.479354</td><td>0.424992</td><td>0.363549</td><td>0.352950</td><td>0.343506</td><td>0.313838</td></tr></tbody></table>

# 4 Numerical study of the stability of solutions of $(GL)_{\infty}$

In [13], Bolley and Helffer proved that the functional $\varepsilon_{\infty}$ is semi-bounded when $h < \frac{1}{\sqrt{2}}$; they have deduced the existence of a solution of the Ginzburg-Landau system corresponding to a global minimum of $\varepsilon_{\infty}$ for any pair $(\kappa, h)$, with $\kappa > 0$ and $0 < h < \frac{1}{\sqrt{2}}$.

The study of the stability of the solutions of (1) - (2) for $h \ge \frac{1}{\sqrt{2}}$ remains an open question. Our investigations, based on the numerical approach of Subsection 3.1 and on a framework involving a generalized eigenvalue problem (see also [18]), will allow us to provide numerical answers to this question.

## 4.1 Generalized eigenvalue problem

Let us start by introducing some notation. For $f_0 \in (0, 1]$, let $(f, A; h)$ be the unique solution of (1) - (3) obtained from Theorem 1.1. Let us set:

$$ M_{f,A} := \begin{pmatrix} 3f^2 + A^2 - 1 & 2fA \\ 2fA & f^2 \end{pmatrix}. \qquad (18) $$

We have studied in [18] the operator $\mathcal{L}_{f,A}$,

$$ \mathcal{L}_{f,A} \begin{pmatrix} g \\ b \end{pmatrix} = \begin{pmatrix} -\kappa^{-2} \frac{d^2}{dx^2} & 0 \\ 0 & -\frac{d^2}{dx^2} \end{pmatrix} \begin{pmatrix} g \\ b \end{pmatrix} + M_{f,A} \begin{pmatrix} g \\ b \end{pmatrix}, \qquad (19) $$
defined on

$$ \mathcal{D} = \{(g, b) \in H^2(\mathbb{R}^+) \times H^2(\mathbb{R}^+) ; g'(0) = 0, b'(0) = 0\}. $$

To check if the pair $(f, A)$ is a local minimum of $\varepsilon_\infty$, we have just to analyze the Hessian of $\varepsilon_\infty$ at $(f, A; h)$. This Hessian is defined on $H^1(\mathbb{R}^+)^2 \times H^1(\mathbb{R}^+)^2$, and it is positive definite if the bottom of the spectrum of the associated operator, $\mathcal{L}_{f,A}$, is strictly positive.

Let us denote by $\sigma(\mathcal{L}_{f,A})$, $\sigma_{ess}(\mathcal{L}_{f,A})$ and $\sigma_{dis}(\mathcal{L}_{f,A}) := \sigma(\mathcal{L}_{f,A}) \setminus \sigma_{ess}(\mathcal{L}_{f,A})$ respectively the spectrum, the essential spectrum and the discrete spectrum of $\mathcal{L}_{f,A}$. The operator $\mathcal{L}_{f,A}$ is self-adjoint and $\sigma_{ess}(\mathcal{L}_{f,A})$ is equal to $[1, +\infty)$. If $f_0 = 1$, the solution of (1) - (3) is $(f, A) = (1, 0)$, and $\sigma_{dis}(\mathcal{L}_{f,A}) = \emptyset$. This particular value of $f_0$ will thus be excluded in our numerical investigations.

Let us denote by $\langle \cdot, \cdot \rangle$ the usual scalar product on $(L^2(\mathbb{R}^+))^2$ and by $\| \cdot \|$ the associated norm. The following statement specifies in particular that, for $f_0 \in (0, 1)$, the discrete spectrum of $\mathcal{L}_{f,A}$ is not empty.

**Proposition 4.1** For $f_0 \in (0, 1)$, let $(f, A; h)$ be the unique solution of $(GL)_\infty$ satisfying (3). For any $\kappa > 0$, there exists $\lambda \in \sigma_{dis}(\mathcal{L}_{f,A})$ such that $\lambda < 1$.

*Proof.* Let us consider the self-adjoint operator $L_1 := \mathcal{L}_{f,A} - \text{id}$ defined on $\mathcal{D}$. We have $\sigma_{ess}(L_1) = [0, +\infty)$. On the other hand, $\sigma_{dis}(L_1) \neq \emptyset$ if and only if there exists $\tilde{\phi} \in \mathcal{D}$ such that $\langle L_1 \tilde{\phi}, \tilde{\phi} \rangle < 0$ (see [32, p. 55]). Let us define $\tilde{\phi} = (g, b)^T := (0, 2e^{-\alpha \kappa x} - e^{-2\alpha \kappa x})^T \in \mathcal{D}$, with $\alpha > 0$ a parameter depending on $f_0 \in (0, 1)$. From this definition, $\int_0^{+\infty} (b'(x))^2 dx = \frac{\alpha \kappa}{3}$ and since we get from (4),

$$ f(x)^2 - 1 \leq -(\cosh(\frac{\kappa x}{\sqrt{2}} + x_0))^{-2} \leq -e^{-2x_0} e^{-\sqrt{2\kappa x}}, \quad \text{where } \tanh(x_0) = f_0, \quad (20) $$

it follows that

$$ \int_{0}^{+\infty} (f(x)^2 - 1)b(x)^2 dx \leq -e^{-2x_0} s(\alpha)\kappa^{-1}, $$

where

$$ s(\alpha) := \frac{22\alpha^2 + 9\sqrt{2}\alpha + 2}{(2\alpha + \sqrt{2})(3\alpha + \sqrt{2})(4\alpha + \sqrt{2})}. $$

In order to get $\tilde{I} := \int_0^{+\infty} ((b'(x))^2 + (f(x)^2 - 1)b(x)^2) dx < 0$, we set $\alpha = \frac{e^{-2x_0}}{10^4}$. Then,

$$ \tilde{I} \leq e^{-2x_0} \left( \frac{\kappa}{3 \times 10^4} - \min_{\alpha \in (0, \frac{1}{10^4}]} s(\alpha) \kappa^{-1} \right), $$

and it is sufficient to choose: $\kappa \leq 10^2$. Thus, as $\langle L_1 \tilde{\phi}, \tilde{\phi} \rangle = \tilde{I}$, it follows that there exists $\kappa_0$, $\kappa_0 = 10^2$, such that, for all $\kappa \leq \kappa_0$, we have $\langle L_1 \tilde{\phi}, \tilde{\phi} \rangle < 0$.

Similarly, in the large $\kappa$-regime, we define $\tilde{\phi} = (g, b)^T := (0, 2e^{-\frac{\beta x}{\kappa^2}} - e^{-\frac{2\beta x}{\kappa^2}})^T$, with $\beta > 0$ a parameter depending on $f_0$. Now, $\int_0^{+\infty} (b'(x))^2 dx = \frac{\beta}{3\kappa^2}$ and from (20),

$$ \int_{0}^{+\infty} (f(x)^2 - 1)b(x)^2 dx \leq -e^{-2x_0} \left( \frac{4}{\frac{2\beta}{\kappa^2} + \sqrt{2}\kappa} - \frac{4}{\frac{3\beta}{\kappa^2} + \sqrt{2}\kappa} + \frac{1}{\frac{4\beta}{\kappa^2} + \sqrt{2}\kappa} \right). $$
By setting $\beta = e^{-2x_0}$, it is sufficient to take $\kappa$ such that: $\frac{1}{3\kappa^2} - \frac{4\beta}{\kappa^2} + \sqrt{2\kappa} < 0$, namely for instance $\kappa \ge 2$, in order to have $\tilde{I} < 0$. It results then that there exists $\tilde{\kappa}_0$, $\tilde{\kappa}_0 = 2$, such that, for $\kappa \ge \tilde{\kappa}_0$, we have $\langle L_1 \tilde{\phi}, \tilde{\phi} \rangle < 0$. As $\tilde{\kappa}_0 \le \kappa_0$, the proof of the statement follows. $\square$

Our approach is here reduced to the numerical inspection of this spectrum, $\sigma_{dis}(\mathcal{L}_{f,A})$, for studying then the stability of the pair $(f, A)$. This approach involves, as in the framework developed in [18], a weak eigenvalue problem associated with $\tilde{\mathcal{L}}_{\tilde{f},\tilde{A}}$ on a truncated domain $[0, d]$, where the homogeneous Dirichlet boundary condition is used at the point $d$. Typically, $d$ takes large values in the physical sense that: $\kappa d \gg 1$. The operator $\tilde{\mathcal{L}}_{\tilde{f},\tilde{A}}$, where $(\tilde{f}, \tilde{A}) := (f_{|[0,d]}, A_{|[0,d]})$,

$$ \tilde{\mathcal{L}}_{\tilde{f},\tilde{A}} \begin{pmatrix} \tilde{g} \\ \tilde{b} \end{pmatrix} = \begin{pmatrix} -\kappa^{-2} \frac{d^2}{dx^2} & 0 \\ 0 & -\frac{d^2}{dx^2} \end{pmatrix} \begin{pmatrix} \tilde{g} \\ \tilde{b} \end{pmatrix} + M_{\tilde{f},\tilde{A}} \begin{pmatrix} \tilde{g} \\ \tilde{b} \end{pmatrix}, $$

is defined on

$$ \tilde{\mathcal{D}} = (H_{00}^2([0, d]))^2; \quad H_{00}^2([0, d]) := \{u \in H^2([0, d]) ; u(d) = 0, u'(0) = 0\}. $$

In the present approach, we shall numerically inspect then the bottom of the discrete spectrum of $\mathcal{L}_{f,A}$ by considering the weak eigenvalue problem associated with $\tilde{\mathcal{L}}_{\tilde{f},\tilde{A}}$, and that consists of finding $(\alpha, (g, b)) \in \mathbb{R} \times \tilde{\mathcal{D}}$ such that: $\forall (u, v) \in \tilde{\mathcal{D}}$,

$$ \left\{ \begin{array}{l} \displaystyle \kappa^{-2} \int_0^d g'(x)u'(x) \, dx + \int_0^d (A(x)^2 + 3f(x)^2 - 1)g(x)u(x) \, dx + 2 \int_0^d A(x)f(x)b(x)u(x) \, dx = \\ \displaystyle \qquad \alpha \int_0^d g(x)u(x) \, dx, \\[2ex] \displaystyle \int_0^d b'(x)v'(x) \, dx + 2 \int_0^d A(x)f(x)g(x)v(x) \, dx + \int_0^d f(x)^2 b(x)v(x) \, dx = \\ \displaystyle \qquad \alpha \int_0^d b(x)v(x) \, dx. \end{array} \right. \qquad (21) $$

Hereafter, the usual scalar product on $(L^2([0, d]))^2$ is represented by $(\cdot, \cdot)_{(L^2([0,d]))^2}$. Also, the (usual) norm associated with any functional space $X$ is denoted by $\| \cdot \|_X$ below.

**Remark 4.1** Let us set $D := \{u \in H^1([0, d]), u(d) = 0\}$, $H := L^2([0, d]) \times L^2([0, d])$, $\tilde{H} := D \times D$; the injection of $\tilde{H}$ in $H$ is compact. Since there exist constants $\gamma > 0$ and $\alpha > 0$ such that $a(U, U) + \gamma \|U\|_H^2 \ge \alpha \|U\|_{\tilde{H}}^2$, $\forall U \in \tilde{H}$, where

$$ a(U, V) := \int_{0}^{d} (\kappa^{-2} u'_{1} v'_{1} + u'_{2} v'_{2} + (A^2 + 3f^2 - 1)u_{1}v_{1} + 2Af(u_{1}v_{2} + u_{2}v_{1}) + f^2 u_{2}v_{2}) dx, $$

the problem that consists of finding a pair $(\lambda, U) \in \mathbb{R} \times \tilde{H}$, $U = (u_1, u_2)^T$, such that:

$$ a(U, V) = \lambda(U, V)_{(L^2([0,d]))^2}, \quad \forall V := (v_1, v_2)^T \in \tilde{H}, $$
admits at least one solution. In fact, from [33, p. 137], there exists $(\lambda_n, U^n) \in \mathbb{R} \times \tilde{H}, U^n = (u_1^n, u_2^n)^T$, with $(\lambda_n)_n$ an increasing sequence, such that $a(U^n, V) = \lambda_n(U^n, V)_{(L^2([0,d]))^2} \ \forall V \in \tilde{H}$. Taking successively $V = (v_1, 0)^T$ and $V = (0, v_2)^T$ in this last variational equation, we get

$$ \left\{ \begin{array}{l} \displaystyle \int_0^d (\kappa^{-2}(u_1^n)'v_1' + (A^2 + 3f^2 - 1)u_1^n v_1 + 2Af u_2^n v_1) \, dx = \lambda_n \int_0^d u_1^n v_1 \, dx \quad \forall v_1 \in D, \\[1.5ex] \displaystyle \int_0^d ((u_2^n)'v_2' + 2Af u_1^n v_2 + f^2 u_2^n v_2) \, dx = \lambda_n \int_0^d u_2^n v_2 \, dx \quad \forall v_2 \in D. \end{array} \right. $$

On the other hand, it results from standard arguments that $U^n \in \tilde{D}$ and satisfies: $\tilde{\mathcal{L}}_{\tilde{f}, \tilde{A}}(U^n) = \lambda_n U^n$.

Following Remark 4.1, the weak problem (21) admits at least one solution.

**Remark 4.2** *From (18) - (19), there exists a real constant C such that: $\langle \mathcal{L}_{f,A}\phi, \phi \rangle \ge C\|\phi\|^2$, $\forall \phi \in (L^2(\mathbb{R}^+))^2$. The operator $\mathcal{L}_{f,A}$ being self-adjoint and its essential spectrum being equal to $[1, +\infty)$, from [32, p. 120], the lowest eigenvalue of $\mathcal{L}_{f,A}$, denoted by $\alpha^\infty(f_0)$, can then be expressed as follows: $\alpha^\infty(f_0) = \inf_{\phi \in \tilde{D}\setminus\{0\}} \frac{\langle \mathcal{L}_{f,A}\phi, \phi \rangle}{\|\phi\|^2}$. Representing similarly the lowest eigenvalue of $\tilde{\mathcal{L}}_{\tilde{f}, \tilde{A}}$ by $\alpha^d(f_0)$, it follows that for any $\epsilon > 0$, there exists $d_0(\epsilon) > 0$ such that for $d \ge d_0(\epsilon)$, we have: $|\alpha^\infty(f_0) - \alpha^d(f_0)| \le \epsilon$.*

As in [18], we also consider here Hermite finite elements (see e.g. [34, Chap. 3]) that appear well suited to an accurate discretization of the space $\tilde{D}$. Let us denote by $\{x_j\}_{1 \le j \le N+1}$ the set of subdivision points of $[0, d]$ where $x_1 = 0, x_{N+1} = d, x_{j+1} = x_j + k_j, N+1$ is the total number of subdivision points of $[0, d]$ and $0 < k_j < 1$. The collection of intervals $I_j := [x_j, x_{j+1}]$, $1 \le j \le N$, will provide a discretization of $[0, d]$. Let us set $k := \max_{1 \le j \le N} \{k_j\}$ and call $\mathcal{T}_k$ such a discretization for which the generic element is an interval denoted $I$. Let us introduce

$$ U_k = \{u_k; u_k, u'_k \in C^0([0, d]), u_k|_I \in P_3, \forall I \in \mathcal{T}_k\}, $$

where $P_3$ is the space of polynomials of degree less than or equal to 3. The discrete space associated with $H_{00}^2([0, d])$ is defined as follows:

$$ U_k^0 = \{u_k \in U_k; u_k(d) = 0, u'_k(0) = 0\}. $$

The discrete problem associated with (21) consists of finding $(\alpha_k, g_k, b_k) \in \mathbb{R} \times U_k^0 \times U_k^0$ such that: $\forall (u_k, v_k) \in U_k^0 \times U_k^0$,

$$ \left\{ \begin{array}{l} \displaystyle \int_0^d g_k'(x) u_k'(x) \, dx + \int_0^d (A(x)^2 + 3f(x)^2 - 1) g_k(x) u_k(x) \, dx + 2 \int_0^d A(x) f(x) b_k(x) u_k(x) \, dx = \\ \displaystyle \qquad \alpha_k \int_0^d g_k(x) u_k(x) \, dx, \\[1.5ex] \displaystyle \int_0^d b_k'(x) v_k'(x) \, dx + 2 \int_0^d A(x) f(x) g_k(x) v_k(x) \, dx + \int_0^d f(x)^2 b_k(x) v_k(x) \, dx = \\ \displaystyle \qquad \alpha_k \int_0^d b_k(x) v_k(x) \, dx. \end{array} \right. \qquad (22) $$
We are thus concerned with (22) in numerical simulations, where A and f are respectively replaced by their numerical approximations $A^{app}$ and $f^{app}$, obtained from the numerical approach of Subsection 3.1. As in [18], we consider the usual Jacobi method (see e.g. [34, Chap. 6]) for the computation of eigenvalues, based on the linear matrix system deriving from (22), where an adaptive Simpson formula is applied for calculating the integral terms.

## 4.2 Numerical results

The numerical study of the stability of the solutions of $(GL)_{\infty}$ is here performed through the numerical inspection of the map $\phi_{k,d}: f_0 \in (0,1) \mapsto \alpha_k(f_0,d)$, where $\alpha_k(f_0,d)$ is the lowest eigenvalue resulting from (22).

We first investigate the weak-$\kappa$ regime. In the experiments associated with the results of Figure 7, we consider $\kappa = 0.1$, and $d = 1000$.

<figure><img src="image_8.png" /><figcaption>Figure 7: Behavior of $\phi_{k,d}$ with respect to values of $f_0$, when $\kappa = 0.1$, $d = 1000$.</figcaption></figure>

For these considered values of $\kappa$ and $d$, there exists a critical value of $f_{0,num}^{eig}$ such that for $f_0 > f_{0,num}^{eig}$, we have $\phi_{k,d}(f_0) > 0$.

We already mention that, in the simulations and following the case, each sequence of values considered for $f_0$ is enriched, between its two consecutive elements providing two eigenvalues of different signs, in order to determine with more accuracy the particular initial value where $\phi_{k,d}$ changes the sign.

Experiments have been performed for examining the influence of the bound $d$ of the truncated domain on the computation of the minimal eigenvalue of the discrete spectrum, in this case where the mesh size $k = 0.05$ is used. By considering $\beta_{d_1,d_2}^k := |\alpha_k(f_0, d_1) - \alpha_k(f_0, d_2)|$, where $d_1, d_2$ are two bounds, we find for instance that the results are quite accurate for $\kappa d = 15$, with thus 150 as the bound of the domain when $\kappa = 0.1$. The curves of Figure 7 are similar to those obtained with such values of $\kappa, d$, and cannot be distinguished from a usual superposition. In fact, for this choice of $\kappa$, it derives from computations that $\beta_{1000,d}^k$ is smaller than $10^{-11}$, when $d = 150, 300, 600, 750$ or 900.
<table><caption>Table 5: Values of $f_{0,num}^{eig}$ in the weak-$\kappa$ regime.</caption><tr><td>κ</td><td>0.1</td><td>0.05</td><td>0.01</td><td>0.004</td><td>0.003</td></tr><tr><td>f<sup>eig</sup><sub>0,num</sub></td><td>0.687382</td><td>0.696754</td><td>0.704931</td><td>0.706161</td><td>0.706376</td></tr></table>

A study of the influence of the mesh size $k$ on the computation of the minimal eigenvalue of the discrete spectrum has also been performed. By considering now $\beta_d^{k_1,k_2} := |\alpha_{k_1}(f_0, d) - \alpha_{k_2}(f_0, d)|$, where $k_1, k_2$ are two mesh sizes, we observe for instance that the results are quite satisfactory with 0.2 as the mesh size, when $d = 150$ and $\kappa = 0.1$; namely, $\beta_d^{0.05,k}$ is smaller than $10^{-4}$, with $k \in \{0.1, 0.2, 0.5\}$. Mainly interested in the change of the sign of $\phi_{k,d}$, we fix $k = 0.2$ for the next computations.

The experiments associated with Figure 8 are performed with $d = 5000$, as smaller values of $\kappa$ are now considered. The results of these experiments, which are not quantitatively influenced when we use a larger value of $d$ in simulations, indicate with respect to $\kappa$ the same remarks as previously. Namely, for each considered value of $\kappa$, it appears that there exists a critical value $f_{0,num}^{eig}$ such that for $f_0 > f_{0,num}^{eig}$, we have $\phi_{k,d}(f_0) > 0$.

<figure><img src="image_9.png" /><figcaption>Figure 8: Behavior of $\phi_{k,d}$ with respect to values of $f_0$, with $d = 5000$ and when $\kappa = 0.1, 0.05, 0.01, 0.003$.</figcaption></figure>

The values obtained from simulations, with respect to $\kappa$, and corresponding to $f_{0,num}^{eig}$, are presented in Table 5.

By comparing these values with those of $f_{0,num}^{sh}$ that were presented in Table 2, we conjecture that there exist $\kappa_0 > 0$ and a map $\theta$, defined on $[0, \kappa_0]$, satisfying $\theta(0) = \frac{1}{\sqrt{2}}$, such that, for any $f_0 > \theta(\kappa)$, the solution $(f, A; h)$ is stable, and for $f_0 \le \theta(\kappa)$, the solution $(f, A; h)$ is unstable.

The numerical results regarding the stability of the solutions of $(GL)_\infty$ are then here in accordance with those obtained in a formal way by Dolgert et al. [20] in the weak-$\kappa$ regime.
We are now concerned with the experiments in the large-$\kappa$ regime. Figure 9 presents the results obtained with $\kappa = 1.3$ and $d = 300$. It appears for these considerations that there exists a critical value $f_{0,\text{num}}^{\text{eig}}$ such that for $f_0 > f_{0,\text{num}}^{\text{eig}}$, we have $\phi_{k,d}(f_0) > 0$.

<figure><img src="image_10.png" /><figcaption>Figure 9: Behavior of $\phi_{k,d}$ with respect to values of $f_0$, when $d = 300$, $\kappa = 1.3$.</figcaption></figure>

A study, in the large-$\kappa$ regime, on the influence of the bound of the truncated domain on the computation of the minimal eigenvalue of the discrete spectrum, has been performed. We observe that the agreement is quite good with $d = 300$; for instance, $\beta_{1000,300}^k$ is smaller than $10^{-7}$, when $\kappa = 1.3$.

We thus fix $d = 300$ for the following experiments. Those associated with Figure 10 are performed with larger values of $\kappa$. It results also, for each of these values, that there exists a critical value $f_{0,\text{num}}^{\text{eig}}$ such that for $f_0 > f_{0,\text{num}}^{\text{eig}}$, we have $\phi_{k,d}(f_0) > 0$.

<figure><img src="image_11.png" /><figcaption>Figure 10: Behavior of $\phi_{k,d}$ with respect to values of $f_0$, when $d = 300$; at left $\kappa = 1.3, 3, 5$, and $\kappa = 9, 11, 15$ at right.</figcaption></figure>

Table 6 reproduces the values obtained from simulations, with respect to $\kappa$, and corresponding to $f_{0,\text{num}}^{\text{eig}}$. A comparison of Tables 4 and 6 shows that these values are close to those of $f_{0,\text{num}}^{\text{sh}}$, in the
<table><caption>Table 6: Values of $f_{0,num}^{eig}$ in the large-$\kappa$ regime.</caption><tr><td>κ</td><td>1.3</td><td>3</td><td>5</td><td>9</td><td>10</td><td>11</td><td>15</td></tr><tr><td>f<sup>eig</sup><sub>0,num</sub></td><td>0.561388</td><td>0.480010</td><td>0.427910</td><td>0.370880</td><td>0.361000</td><td>0.351080</td><td>0.324300</td></tr></table>

<table><caption>Table 7: Values of $f_{0,num}^{sh}$ and $f_{0,num}^{eig}$ in the intermediary-$\kappa$ regime.</caption><tr><td>κ</td><td>0.6</td><td>0.7</td><td>0.707</td><td>0.8</td><td>0.9</td></tr><tr><td>f<sup>sh</sup><sub>0,num</sub></td><td>0.620053</td><td>0.609918</td><td>0.609242</td><td>0.600493</td><td>0.591678</td></tr><tr><td>f<sup>eig</sup><sub>0,num</sub></td><td>0.620054</td><td>0.609928</td><td>0.609247</td><td>0.600512</td><td>0.591713</td></tr></table>

large-$\kappa$ regime.

Since the same remark was observed in the weak-$\kappa$ regime, we are then led to conjecture that: in each regime of $\kappa$, $\sigma_{\kappa}(f_0) = h^{sh}(\kappa)$ if and only if $\phi_{k,d}(f_0) = 0$.

The results presented in Table 7 suggest that this remark appears also to be satisfied in the context of intermediary values of $\kappa$.

The values of $f_{0,num}^{eig}$ become, very slowly, smaller when $\kappa$ increases; the consideration of $\kappa = 21$ for example provides $f_{0,num}^{eig} = 0.297935$ which is not significantly smaller than the value that was obtained for $\kappa = 15$. Our results suggest in the large-$\kappa$ regime that there exist $\tilde{\kappa}_0 > 0$ and a map $\tilde{\theta}$ decreasing on $[\tilde{\kappa}_0, +\infty)$, such that for all $f_0 \le \tilde{\theta}(\kappa)$ the solution $(f, A; h)$ is unstable, and for all $f_0 > \tilde{\theta}(\kappa)$ the solution $(f, A; h)$ is stable.

In each regime of $\kappa$, these results suggest also that $\phi_{k,d}$ is continuous on $(0, 1)$, and admits a unique minimum on $(0, 1)$ achieved at $\gamma(\kappa)$. Moreover, we can conjecture that $\phi_{k,d}$ is decreasing on $(0, \gamma(\kappa)]$, and increasing on $[\gamma(\kappa), 1)$.

# 5 Conclusions

A numerical inspection of certain maps and conjectures associated with the solutions of the half-space Ginzburg-Landau system has been performed in the weak and large-$\kappa$ regimes. These conjectures involve the superheating field $h^{sh}$; namely, the asymptotic expansion of $h^{sh}$ proposed by Chapman [30], in the large-$\kappa$ regime, as well as the estimate (cf. [22]) involving $h^{sh}$ in all the regimes of $\kappa$ appear numerically satisfied. On the other hand, we have numerically studied the stability of solutions of this system in all the regimes of $\kappa$. In the weak-$\kappa$ regime, our numerical results suggest that there exist $\kappa_0 > 0$ and a map $\theta$ defined on $[0, \kappa_0]$, with $\theta(\kappa) \in (0, 1)$, $\theta(0) = \frac{1}{\sqrt{2}}$, such that for all $f_0 \le \theta(\kappa)$ the solution $(f, A; h)$ is unstable, and for all $f_0 > \theta(\kappa)$ the solution $(f, A; h)$ is stable. Typically, these results extend those that we had obtained in [18] from the framework based on a localization of solutions of the system $(GL)_{\infty}$ satisfying the de Gennes condition. In the large-$\kappa$ regime, our numerical results suggest that there exist $\tilde{\kappa}_0 > 0$ and a map $\tilde{\theta}$ decreasing on $[\tilde{\kappa}_0, +\infty)$, with $\tilde{\theta}(\kappa) \in (0, 1)$, such that for all $f_0 \le \tilde{\theta}(\kappa)$ the solution $(f, A; h)$ is unstable, and for all $f_0 > \tilde{\theta}(\kappa)$ the solution $(f, A; h)$ is stable.
the solution $(f, A; h)$ is stable.

Our numerical investigations lead us to conjecture that $\sigma_{\kappa}(f_{0,c}) = h^{sh}(\kappa)$ if and only if 0 is the lowest eigenvalue of $\mathcal{L}_{f,A}$, where $f(0) = f_{0,c}$. In other words, such a value $f_{0,c}$ corresponds to the initial datum for which the change of stability occurs.

The results obtained here provide the numerical answers to several questions that remain open from a theoretical point of view and concern therefore the stability of solutions of the half-space Ginzburg-Landau system as well as the study of the two conjectures involving (in the large-$\kappa$ regime and in all the regimes of $\kappa$ respectively) the field $h^{sh}$. The theoretical proof or a numerical inspection of the asymptotic expansion of $h^{sh}$ in powers of $\kappa^{1/2}$, introduced in [29] formally, at any order, and extending thus in the weak-$\kappa$ regime the Parr formula (see (7)), remains also to be performed. Another numerical work is the one that would investigate, in the large-$\kappa$ regime, the extension of the asymptotic expansion of $h^{sh}$ proposed by Chapman (see (17)), by describing then the associated third term.

**Acknowledgments.** The authors thank Professor C. Bolley for her useful comments.

## References

1. Bolley C, Helffer B (2002) Global superheating field for superconductors in a large bounded interval. Physica D 172:162-189

2. Bolley C, Helffer B (1996) Rigorous results for the Ginzburg-Landau equations associated with a superconducting film in the weak $\kappa$-limit. Reviews in Math Physics 8:43-83

3. Ginzburg VL (1955) On the theory of superconductivity. Il Nuovo Cimento 2(6):1234-1250

4. Ginzburg VL (1958) On the destruction and the onset of superconductivity in a magnetic field. Zh Eksp Teor Fiz 34:113-125

5. Ginzburg VL, Landau LD (1950) On the theory of superconductivity. Zh Eksp Teor Fiz 20:1064-1082

6. Pretzl K (1993) Superconducting granule detectors. Journal of Low Temperature Physics 93:439-448

7. Tinkham M (1975) Introduction to superconductivity. McGraw-Hill, New York

8. Chapman SJ, Howison SD, Ockendon JR (1992) Macroscopic models for superconductivity. SIAM Rev 34:529-560

9. Kwong MK (1995) On the one-dimensional Ginzburg-Landau BVP's. Differential Integral Equations 8:1395-1405
10. Yang Y (1990) The Ginzburg-Landau equations for a superconductivity film and the Meissner effect. J Math Phys 31:1284-1289

11. Aftalion A, Troy WC (1999) On the solutions of the one-dimensional Ginzburg-Landau equations for superconductivity. Physica D 132:214-232

12. Yang Y (1990) Boundary value problems of the Ginzburg-Landau equation. Proc Roy Soc Edin Sect A 114:355-365

13. Bolley C, Helffer B (1997) The Ginzburg-Landau equations in a semi-infinite superconducting film in the large $\kappa$ limit. Eur J Appl Math 8:347-367

14. Parks RD (ed) (1969) Superconductivity: Volume 2. Marcel Dekker, New York

15. Bolley C, Foucher F, Helffer B (2000) Superheating field for the Ginzburg-Landau equations in the case of a large bounded interval. J Math Phys 41:7263-7289

16. Bolley C, Del Castillo P (2001) Existence and uniqueness for the half-space Ginzburg-Landau model. Nonlinear Analysis 47:135-146

17. Bolley C, Helffer B (2006) On the convergence of an alternative shooting method on the half-space Ginzburg-Landau system. M³AS 16:1527-1558

18. Del Castillo P, Mefire SM (2005) Some aspects of the one-dimensional Ginzburg-Landau system in the weak-$\kappa$ limit and numerical study of the stability of solutions. Far East J Appl Math 19:81-123

19. De Gennes PG (1966) Superconductivity. Selected topics in solid state physics and theoretical physics. In: Proceedings of the 8th Latin American School of Physics, Caracas, 1966

20. Dolgert AJ, Di Bartolo SJ, Dorsey AT (1996) Superheating fields of superconductors: Asymptotic analysis and numerical results. Physical Review B 53:5650-5660. Erratum: (1997) Superheating fields of superconductors: Asymptotic analysis and numerical results. Physical Review B 53:5650-5660. Physical Review B 56:2883-2883

21. Bolley C, Helffer B (1997) Proof of the de Gennes formula for the superheating field in the weak-$\kappa$ limit. Ann de l'Inst H Poincaré 14:597-613

22. Bolley C, Helffer B (2002) A priori estimates for Ginzburg-Landau solutions. In: Berestycki H, Pomeau Y (eds) Nonlinear PDE's in Condensed Matter and Reactive Flows. Conference in Institut Scientifique de Cargèse, 1999. NATO Science Series, vol 569. Kluwer, pp355-373

23. Del Castillo P (2005) The Parr formula for the superheating field in a semi-infinite film. J Math Phys 46(5):053513-1-053513-22. Publisher's note: (2005) The Parr formula for the superheating field in a semi-infinite film [J Math Phys 46(5):053513-1-053513-22] by P. Del Castillo. J Math Phys 46(7):079902-1-079902-1

24. Del Castillo P (2002) Two terms in the lower bound for the superheating field in a semi-infinite film in the weak-$\kappa$ limit. Eur J Appl Math 13:519-547
25. Bolley C, Helffer B (1997) Superheating in a semi-infinite film in the weak-$\kappa$ limit: numerical results and approximate models. Math Model Num Anal 31:121-165

26. Crouzeix M, Mignot AL (1984) Analyse numérique des équations différentielles. Masson, Paris

27. Quarteroni A, Sacco R, Saleri F (2000) Numerical mathematics. Springer Verlag, New York

28. Parr H (1976) Superconductive superheating field for finite $\kappa$. Z Physik B 25:359-361

29. Del Castillo P (2002) Expansion for the superheating field in a semi-infinite film in the weak-$\kappa$ limit. Math Model Num Anal 36:971-993.

30. Chapman SJ (1995) Superheating field of type II superconductors. SIAM J Appl Math 55:1233-1258

31. Helffer B, Weissler FB (1998) On a family of solutions of the second Painlevé equation related to superconductivity. Eur J Appl Math 9:223-243

32. Hislop PD, Sigal IM (1996) Introduction to spectral theory: with applications to Schrödinger operators. Springer Verlag, New York

33. Raviart PA, Thomas JM (1988) Introduction à l'analyse numérique des équations aux dérivées partielles. Masson, Paris

34. Ciarlet PG (1998) Introduction à l'analyse numérique matricielle et à l'optimisation. Dunod, Paris