Sweeping Algorithms for Inverting the Discrete Ginzburg-Landau Operator^{∗}

Man Kam Kwong

Mathematics and Computer Science Division Argonne National Laboratory Argonne, IL 60439-4844

Abstract

The Ginzburg-Landau equations we study arise in the modeling of superconductivity. One widely used method of discretizing the equations together with the associated periodic boundary conditions, in the case of a rectangle of dimension 2, leads to a five-point stencil. Solving the system means inverting a sparse matrix of dimension N^{2}, where N is the number of grid points on each side of the rectangle. We propose a method that is similar to the shooting technique in the numerical solution of ordinary differential equations. For small N, the method requires inverting a full matrix of dimension 2N. When N is large, an iterative procedure combining partial sweeping and the technique of divide and conquer (domain decomposition) is appropriate.

AMS(MOS) Subject Classification. Primary 65F05, 65N22. Secondary 82D55.

Key Words and Phrases. Numerical computation, boundary value problem, Ginzburg-Landau equations, superconductivity.

Proposed Running Head. Inverting the Ginzburg-Landau operator

∗This work was supported by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S. Department of Energy, under Contract W-31-109-Eng-38.

1 Introduction

The term Ginzburg-Landau equations has been applied to several sets of differential equations arising in diversified applications. The mathematical properties of the various sets, however, can differ widely. The particular equations we are interested in here were proposed by Ginzburg and Landau [13] in 1950 to model the phenomenon of superconductivity. They are derived as the first variational equations for the minimization of a certain free energy functional. One of the important goals of research is the understanding of the behavior of solutions satisfying given boundary conditions. The boundary conditions we shall use are the so-called quasi-periodic conditions.

Because the solutions do not have simple analytical representations, numerical computation becomes an important tool. A nonstandard finitedifference discretization has been used by researchers in lattice gauge field theory, in order to preserve a discrete version of gauge invariance. A brief discussion of the Ginzburg-Landau theory and the discretization method will be given in Section 2.

Numerical methods that have been tried on two-dimensional cases include simulated annealing [8] and relaxation [19] [16] (which is equivalent to solving a time-dependent version of the Ginzburg-Landau equations using a forward Euler scheme of integration). Although these methods are easy to understand and to implement, they are not efficient and require an enormous amount of computer time. Recently, Wright and Plassmann succeeded in using advanced optimization and numerical linear algebra techniques to efficiently compute the solutions. In collaboration with materials scientists Garner, Spanbauer, Benedek, and Strandburg [12], they were able to extend their algorithm to tackle a three-dimensional version that models multilayered superconductors. A different approach using finite element methods is adopted by Du, Maxburger, and Peterson [9] [10].

The more sophisticated methods require the solution of a large system of nonlinear algebraic equations, which is accomplished with an iterative Newton method. Within each Newton step, however, there is still a large system of linear equations to be solved. In the work reported in [12], the linear system is treated with yet another iterative scheme. The purpose of this paper is to propose an alternative method for solving an important subsystem of the linear system.

Two physical quantities, both functions of position, are involved in the modeling: a complex scalar order parameter ψ and a real-valued vector potential A. The Ginzburg-Landau equations are, respectively, differential equations governing the spatial variation of ψ and A. From the point of view of numerical computation, the equation for ψ is harder to solve, and that is our main concern in this paper.

We consider the two-dimensional case and assume that the part of the superconductor under investigation is represented by a rectangle. As usual, the rectangle is discretized into N^{2} lattice points. Whether standard or nonstandard finite difference is used, the discretized equation for ψ consists of a linear part (strictly speaking, it is linear in ψ, but not in A), which is a fivepoint stencil on the lattice points, and a relatively simple nonlinear function of ψ. The linear part, which we shall call the Ginzburg-Landau operator, can be considered as a modification of the classical Laplace operator.

If one assumes A to be given, for instance, from the most recent iterative step, the equation can be solved for ψ by using a Newton method, provided that one can invert the Ginzburg-Landau operator. In theory, the operator corresponds to a sparse matrix acting on the set of N^{2} unknown ψ values at all the lattice points. However, direct methods for inverting an N^{2} × N^{2} matrix, even though it is sparse, are not practical for large N. Plassmann and Wright employ an iterative algorithm (nested within each Newton step) for this purpose. In theory, our approach inverts the Ginzburg-Landau operator directly, using an idea borrowed from the shooting method for solving ordinary differential equations; see, for example, [3]. We reduce the problem to the inversion of a (full) 2N × 2N matrix. The inverse Ginzburg-Landau image of any given matrix can then be found without the explicit construction of the inverse of the Ginzburg-Landau operator. For small N, the method works beautifully. An obvious advantage of this method is that the inverse of the reduced matrix can be retained for use in subsequent Newton steps, whereas in using an iterative method, one has to start from scratch every time. When N is large, however, the 2N × 2N matrix involved is extremely illconditioned, and conventional methods of inversion are inadequate. In Section 4, we discuss how the idea of divide and conquer can be incorporated. The technique is also the same as that of domain decomposition. The resulting algorithm is particularly suited for parallel computation. We also discuss how large errors can arise at boundaries between adjoining subdomains and what can be done to smooth them out.

A three-dimensional model of a superconductor requires a three dimensional lattice with N^{3} points and a seven-point stencil. The sweeping method can also be extended to reduce the inversion of the stencil to that of a full 2N^{2} × 2N^{2} matrix, instead of the inversion of a sparse N^{3} × N^{3} matrix. How the reduction can be used effectively will be the object of a future investigation.

2 The Ginzburg-Landau Model of Superconductivity

Many good monographs exist on the theory of superconductivity [7] [17] [18]. Several recent survey articles [5] [9] [10] [15] have also been written with mathematicians in mind; these should be consulted for further details on the physics, experimental observations, and alternative models. We describe here only the necessary concepts to make our mathematical problem understandable.

Superconductivity was discovered in 1911 by Kammerlingh Onnes, who observed experimentally that certain metals could lose their electric resistance completely when cooled below a critical temperature T_{c}, characteristic of the metal. He also observed that, while at a temperature below T_{c}, a sufficiently large external magnetic field could revert the metal from the superconducting state to the normal state.

Ginzburg and Landau in 1950 [13] explained the phenomenon as a new type of phase transition, similar to the familiar solid-liquid or liquid-gas phase change of matter. In the case of superconductivity, no outward change in appearance is discernible, but fundamental changes do occur at the atomic level, in particular, from the existence to the nonexistence of superconducting electrons, the so-called Cooper pairs. By drawing an analogy with the theory of classical phase transitions, Ginzburg and Landau obtained their famous model. It is hypothesized that two quantities, ψ and A, which represent the physical state of the metal, must minimize a certain functional, the Helmholtz free energy. A nonlinear term in the functional depends on the temperature T. At T > T_{c}, the normal state for the metal is the minimizer for the functional. Below T_{c}, the normal state is only a stationary point. A “superconducting” minimizer appears at some “distance” away from the normal state.

The quantity ψ is a complex scalar function of position and is usually referred to as the order parameter. It plays a role similar to the wave function in quantum mechanics; the square of its absolute value, |ψ|^{2}, represents the local density of superconducting electrons. The second quantity is the realvalued three-dimensional vector magnetic potential A used in the classical theory of electromagnetism. The value of the magnetic field is given by the Maxwell equation as h = ∇ × A, and the supercurrent is a multiple of ∇×h = ∇×∇×A. Various other observable physical quantities such as the specific heat can likewise be calculated once ψ and A are known. In other words, the electromagnetic state of the metal is determined by ψ and A. We shall see below, however, that the converse is not true.

In the simplest situation, the Helmholtz energy consists of three parts: the condensation energy F (ψ), which is a nonlinear function of the order parameter and temperature; the kinetic energy, which represents the interactions between the vector potential and the order parameter; and the field energy, which depends on h = ∇ × A alone. The form of F (ψ), modulo an additive constant, was postulated by Ginzburg and Landau to be

F (ψ) = α(T)|ψ|^{2} + \frac{1}{2}β(T)|ψ|^{4},  (2.1)

where the coefficients α(T) and β(T) are constant throughout a homogeneous superconductor once T is fixed. When T > T_{c}, α(T) > 0, and when T < T_{c}, α(T) < 0, while for all T, β(T) > 0. The other two energy terms are standard in classical field theory.

The Helmholtz energy can also be written for a general inhomogeneous and anisotropic material, which is the case for the most recently discovered high-T_{c} ceramic materials. Some of these are made up of interlacing layers of superconducting and nonsuperconducting materials, forming a stack of Josephson junctions. The optimization procedure remains the same for such materials as for the simpler homogeneous isotropic metal. The difference is in the numbers of variables and in the complexity of the various coefficients involved. For the purpose of simple elucidation of the computation principle, we restrict ourselves to the simpler case.

Even for a homogeneous isotropic superconductor, the original form of the energy functional involves various physical constants such as the Planck’s constant, the charge and mass of the electron, and the constant π. As far as the mathematics is concerned, most of these constants can be hidden by choosing new scales for the various axes and by scaling the quantities ψ and A, to result in the so-called dimensionless energy functional

G(ψ, A) = ∫_{Ω} (−|ψ|^{2} + \frac{1}{2}|ψ|^{4} + |(∇ − iA) ψ|^{2} + κ^{2} |∇ × A|^{2}) dΩ, (2.2)

where i = \sqrt{−1}.

We use the same symbols ψ and A in the above formula for the functional, but one must bear in mind that the original quantities that correspond to real physical dimensions are to be recovered after a proper scaling. Now |ψ(x)|^{2} signifies the degree of superconducting activity at the location x. In particular, |ψ(x)| = 0 means lack of superconducting current and the point x is said be in the normal state, while |ψ(x)| = 1 signifies perfect superconductivity.

The temperature T no longer appears explicitly in the condensation energy portion of the functional since the coefficients α and β have been reduced to unity. However, the effect of T on α and β is manifested through the choice of the scaling size, which in turn affects the coefficient κ = κ(T) in the last term.

By using a different choice of scaling, one can make κ appear in the kinetic term (replacing (∇+iA) by (∇/κ+iA)), instead of in the last term. The latter is indeed a more commonly used form of the functional; see, for example, [9]. The form (2.2) has been used in [8]. For our purpose, (2.2) is more convenient since the constant κ will not appear in the variational equation governing ψ, which is the main object of our discussion in this paper. Instead, κ appears in the differential equation governing A; it affects ψ indirectly through the coupling of the equations.

For fixed T, the “constant” κ is characteristic of the superconducting material. There is a critical value of κ, namely, κ_{c} = 1/\sqrt{2}. When κ < 1/\sqrt{2} for all T, the superconductor is said to be of Type I; all of the superconductors discovered before 1965 fall into this category. An important property of this class of superconductor is that the magnetic field is expelled from the interior of the material (the Meissner-Ochsenfeld paramagnetic property) and that the supercurrent is confined to a very thin layer at the boundary. This means that both ψ and A decay rapidly to 0 towards the interior of the superconductor. The solutions obtained in solving the Ginzburg-Landau model confirm this phenomenon, leading to wide acceptance of the model. Such Type I superconductors are of limited practical value, however.

In 1957, Abrikosov made a theoretical study of the Ginzburg-Landau equations for κ > 1/\sqrt{2} and discovered that the solutions can behave differently. When the external magnetic field is within an appropriate range (between the so-called H_{c1} and H_{c2} critical fields), the solutions are no longer small in the interior of the superconductor. Instead, |ψ| varies in some periodic pattern (boundary and corner effects cause slight distortions to the periodicity). Furthermore, ψ vanishes at points that form a lattice of a regular pattern. At each of these points, there is no superconducting activity, but in a neighborhood of each, superconducting electrons are moving in circular paths around the point much like vortices in fluid dynamics. By analogy, each such point is said to be the center of a vortex. The material is said to be in a mixed state, since superconducting currents occur almost everywhere in the region, the exceptions being at the vortex centers where the metal is (at least locally) in a normal state. Abrikosov’s prediction was confirmed about ten years later when such Type II superconductors were discovered experimentally. The ability of the magnetic field to penetrate into the interior of Type II materials opens up the possibility of practical applications of superconductors. The recent discovery of high-T_{c} Type II superconductors makes this goal even more plausible, and, as a result, the study of the Ginzburg-Landau equations has surged.

One must bear in mind, however, that the Ginzburg-Landau theory was first proposed as a phenomenological model, implying that its ultimate validity depends on the agreement of theoretical and experimental results. Even though in 1959 Gor’kov [14] was able to derive the Ginzburg-Landau model from the more fundamental BCS (Bardeen, Cooper, and Schrieffer [2]) microscopic theory, some idealization and approximations are adopted in the derivation. For example, the temperature is assumed to be close to T_{c}, and the spatial variations of ψ and A are slow. Furthermore, more recent theories suggest that other physical entities, such as random thermal noise and impurity sites that occur randomly throughout the material, play some roles in the complete description of the behavior of high-T_{c} superconductors. If these theories are correct, the original Ginzburg-Landau model would have to be modified by including such physical entities. In the extreme situation, the model might have to be replaced by a completely new one. In any case, a critical study of the Ginzburg-Landau equations is necessary to determine which action is to be taken.

The mathematical problem is not yet well defined without specifying boundary conditions for ψ and A. If the domain Ω represents the entire piece of superconductor, the physical boundary conditions are given by the facts that no current flows from the interior of the superconductor across the boundary into the ambient space and that the magnetic field at boundary points coincides with the external field. These boundary conditions have been used in various simulations; see, for example, [16] [11]. However, such attempts to simulate an entire piece of a superconductor are not realistic. Not only is the resolution coarse, but the size of the domains is unrealistically small in order to limit the number of vortices. The computational task providing sufficient resolution to include all the vortices is beyond the capacity of the the largest computer available.

Experimental findings revealed a periodic pattern of vortices at some distance away from the boundary. The vortices form a regular triangular (some call it hexagonal) lattice (if one ignores some slight aberrations that can be explained by the presence of impurities). A different approach to compute the Abrikosov vortices numerically is to simulate just an interior portion of the superconductor. The domain Ω is chosen to give a core pattern that can reproduce the periodic vortex pattern by tessellation: The boundary conditions on ∂Ω are then of periodic type. Periodicity is imposed, however, on measurable physical quantities, such as current and magnetic field, but not on ψ and A. This is a consequence of the fact, already alluded to, that the order parameter and the vector potential are not uniquely determined by the physical state. This is an intrinsic characteristic of gauge field theory.

Given ψ, and A, one can pick any real-valued function of position χ(x), called a gauge, and form a new pair

\overline{ψ} = ψe^{iχ(x)},  (2.3)
\overline{A} = A + ∇χ(x).  (2.4)

This mapping is called a gauge transformation. One can verify that the new pair gives the same superconducting current density, the same magnetic field, and the same current as the original pair. Furthermore, the Helmholtz energy computed by using both pairs is the same:

G(\overline{ψ}, \overline{A}) = G(ψ, A).  (2.5)

The energy functional is said to be gauge invariant. The two pairs of potentials are considered to be equivalent representations of the same electromagnetic state of the material. By varying the gauge, one obtains an entire class of equivalent representations.

As a result, the requirement that measurable physical quantities be identical at opposite boundaries points of Ω does not mean that ψ and A must agree at the same points, but only that ψ and A at one boundary point be gauge equivalent to ψ and A at the other boundary point. A thorough explanation of the implication of this statement can be found in [10].

Even after taking the gauge invariance into consideration, there is still some arbitrariness in choosing the boundary conditions. We follow Doria, Gubernatis, and Rainer [8] in using the so-called quasi-periodic boundary conditions.

To be more specific, we take Ω to be the rectangle [0, L_{x}] × [0, L_{y}] which contains n vortices, and assume that A = (A, B, 0) has essentially only two components. The number of vortices depends on how large a core domain we would like to take. We can also pack more vortices into a domain of a given size by increasing the external magnetic field, namely, by increasing the value of g in the following formulas. The boundary conditions are then expressed as

ψ(0, y) = ψ(L_{x}, y),  (2.6)
ψ(x, 0) = ψ(x, L_{y})e^{igx},  (2.7)
A(0, y) = A(L_{x}, y),  (2.8)
A(x, 0) = A(x, L_{y}) + g,  (2.9)
B(0, y) = B(L_{x}, y),  (2.10)
B(x, 0) = B(x, L_{y}),  (2.11)

where

g = \frac{2nπ}{L_{x}} ,  (2.12)

together with the requirement that the first derivatives of the functions involved be continuous after wrapping around across the boundary. Note that our conditions are different from those in [8], but only because we have made a 90 degree rotation in orientation. We find that the computer code is easier to write with this modification.

Du, Gunzburger, and Peterson in [10] used as the core domain, Ω, a parallelogram formed by the centers of four vortices. The quasi-periodic boundary conditions are then imposed on opposite sides of the parallelogram. We shall see later that we can still treat their model by using a rectangular domain.

We next discretize the problem for the purpose of computation. We use a uniform grid of N × N points, not including the points on the top and right-hand boundaries of the rectangle. The grid spacings are then

h_{x} = L_{x}/N and h_{y} = L_{y}/N  (2.13)

in the x and y directions, respectively. If the derivatives are approximated by forward differences, the discretized Helmholtz functional takes the form

G_{d} = \frac{L_{x}L_{y}}{N^{2}} ∑ (−|ψ|^{2} + \frac{1}{2}|ψ|^{4} + |\frac{ψ^{→} − ψ}{h_{x}} + iAψ|^{2}+ |\frac{ψ^{↑} − ψ}{h_{y}} + iBψ|^{2}+

+ κ^{2} |\frac{A^{↑} − A}{h_{y}} − \frac{B^{→} − B}{h_{x}} |^{2}) ,  (2.14)

where we use the arrow notation to denote the value of a function at an appropriate neighboring point. Using central differences, we can improve the accuracy of the approximation, but risk introducing false solutions. Another approach is to discretize not the functional itself but the Ginzburg-Landau equations obtained from taking the first variation of the functional. We do not pursue this alternative but adopt a nonstandard discretization.

A major disadvantage of using (2.14) is that the discretized form of the gauge invariance is only approximately satisfied. Two gauge-equivalent pairs do not necessarily yield identical values for G_{d}. To remedy this, gauge field theorists have long been using a different approximation (see [4]) for the kinetic energy term:

G = \frac{L_{x}L_{y}}{N^{2}} ∑ (−|ψ|^{2} + \frac{1}{2}|ψ|^{4} + |\frac{ψ^{→} − e^{iAh_{x}}ψ}{h_{x}} |^{2} + |\frac{ψ^{↑} − e^{iBh_{y}}ψ}{h_{y}} |^{2} +

+ κ^{2} |\frac{A^{↑} − A}{h_{y}} − \frac{B^{→} − B}{h_{x}} |^{2}) .  (2.15)

This functional is invariant under the discretized gauge transform

\overline{ψ} = ψe^{iχ(x)},  (2.16)
\overline{A} = A + \frac{χ^{→} − χ}{h_{x}},  (2.17)

\overline{B} = B + \frac{χ^{↑} − χ}{h_{y}}.  (2.18)

Taking Taylor expansions, we see that G and G_{d} agree up to first-order terms, as max |A|h_{x}, |B|h_{y} → 0. Another advantage of using (2.15) is that the discrete Ginzburg-Landau equations it leads to are much simpler than that yielded by (2.14). As long as we keep both |A|h_{x} and |B|h_{y} sufficiently small, we may as well choose G over G_{d}.

Another device used by gauge theorists appears at first to be hard to justify. The discretized vector potential A (equivalently, the functions A and B), has been called a bond (or link) variable and is defined not at the grid points, but on the bonds (the short line segments) connecting the grid points. More precisely, the value A lives on the bond between p and p^{→} and B lives on the bond between p and p^{↑}. Rather than following this approach and trying to define precisely what living on the bond means, we find that theoretically it makes no difference if we assume that A and B are defined at the lattice points. We shall see below, however, that associating A and B with the appropriate bonds is a useful way to remember the correct coefficients associated with the five-point stencil representing the discrete Ginzburg-Landau operator.

Our minimization functional depends, in fact, on four real functions, the real and imaginary parts of ψ, A, and B. The discrete Ginzburg-Landau equations are obtained by setting the first variation with respect to each of these to zero. Straightforward computation shows that the first two equations can be combined into one on ψ:

\frac{e^{iA^{←}h_{x}}ψ^{←} − 2ψ + e^{−iAh_{x}}ψ^{→}}{h_{x}^{2}} + \frac{e^{iB^{↓}h_{y}}ψ^{↓} − 2ψ + e^{−iBh_{y}}ψ^{↑}}{h_{y}^{2}} +

+ (1 − |ψ|^{2}) ψ = 0.  (2.19)

The other two equations are simpler to solve from the point of view of numerical computation and will not be studied in this paper.

The first two terms in (2.19) are a five-point stencil: a linear combination of ψ at a point with those at its four neighbors (left, right, up, and down). This is analogous to the standard finite-difference discretization of the Laplacian Δψ, except that the coefficients involved are not the familiar {1, −2, 1} but some variable values depending on A and B. When it comes to boundary points, the five-point stencil in (2.19) has to be interpreted in light of the boundary conditions. For a point on the left boundary of Ω, its left neighbor becomes the point on the right boundary on the same row, and A^{←} becomes A(:, N). There is a wraparound action involved. For a point on the upper boundary of Ω, its upper neighbor is the point on the lower boundary in the same vertical column. However, the coefficient for the upper neighbor becomes e^{iBh_{y}}e^{igx}. Similarly, the coefficient for the “lower” neighbor of a point on the lower boundary is e^{iB(N,:)h_{y}}e^{−igx}.

The Ginzburg-Landau equations are only necessary conditions for a pair (ψ, A) to minimize G. A solution to the equations may not be a global minimizer, however. It can be a local minimizer or even an unstable saddle point. Nevertheless, if one is lucky enough to start with a guess close to the minimizer, a Newton method can yield a sequence of pairs converging to the desired minimizer. Each Newton step requires solving the equation

Lψ_{n} = −(1 − |ψ_{n−1}|^{2})ψ_{n−1}  (2.20)

for ψ_{n}, together with the two other equations in A and B. We have used L to denote the discrete Ginzburg-Landau operator, and ψ_{n−1} the computed value of ψ from the last iterative step.

The approach used by Garner et al. in [12] is based on a conjugate gradient technique in optimization theory, but the nature of each iterative step is similar to that described above and involves the inversion of a linear system that encompasses the discrete Ginzburg-Landau operator.

An efficient method for solving (2.20) (equivalently, for inverting L) can therefore result in a faster scheme for solving the Ginzburg-Landau model. The sweeping method we present in the next Section is devised with this objective in mind.

Our method is equally applicable to operators obtained from discretizing other models.

If we start with (2.14) instead of (2.15), we arrive at a different set of discrete Ginzburg-Landau equations. The one corresponding to (2.19) is

\frac{ψ^{←} − 2ψ + ψ^{→}}{h_{x}^{2}} + \frac{ψ^{↓} − 2ψ + ψ^{↑}}{h_{y}^{2}} − (A^{2} + B^{2})ψ

+\frac{i(A^{←}ψ^{←} − Aψ^{→})}{h_{x}} + \frac{i(B^{↓}ψ^{↓} − Bψ^{↑})}{h_{y}}

+ (1 − |ψ|^{2}) ψ = 0.  (2.21)

The equation looks more complicated, but the linear part is still a five-point stencil, connecting the current point with its four neighbors. Our sweeping method therefore works.

We next look at the oblique model used by Du, Gunzburger, and Peterson [10]. The domain Ω has a horizontal edge and an adjacent edge at an angle of 60^{∘} with the former. Du et al. used the finite element approach which aligns more closely with standard finite difference discretization. We can apply the nonstandard finite difference method to the oblique model. Instead of using a parallelogram grid to mimic the shape of the domain Ω, we find it more convenient to reformulate the problem on a rectangular grid. We cut off a triangle from the right-hand side of the parallelogram, move it to the left to form a rectangle, and then discretize it in the usual way. We choose the grid spacings in the x and y directions appropriately so that the oblique edge of Ω corresponds to an oblique line through appropriate grid points. Interior points inside the grid as well as those on the left and right boundaries satisfy the same Ginzburg-Landau five-point stencil given by (2.19). Those on the upper and lower boundaries, however, satisfy a similar stencil with upper and lower neighbors replaced by points lying on the opposite boundary determined by drawing a line parallel to the oblique edge of Ω. For this type of stencil, the upward/downward version of our sweeping algorithm is applicable.

3 The Sweeping Method

In this section, we describe the theory of the sweeping method. The method is suitable for a general five-point stencil in which the coefficient for the “right” neighbor is never zero. A modified algorithm works when “right” is replaced by “left,” “up,” or “down,” In a simpler form, it can be used to invert banded, almost triangular, or cyclic matrices. Let us start by using a cyclic tridiagonal system as an example. Such systems arise in the discretization of periodic boundary value problems of second-order differential equations.

We would like to solve for the unknown vector x in the equation

Kx = ( \begin{matrix}−a_{1} & 1 & & & −b_{1} \\ −b_{2} & −a_{2} & 1 \\ & . & . & . \\ & & . & . & 1 \\ 1 & & & −b_{n} & −a_{n} \end{matrix} ) ( \begin{matrix}x_{1} \\ x_{2} \\ . \\ . \\ x_{n} \end{matrix} ) = ( \begin{matrix}c_{1} \\ c_{2} \\ . \\ . \\ c_{n} \end{matrix} ) = c. (3.1)

We assume that the matrix K is nonsingular. The method works equally well even if the entries in K below the diagonal are nonzero.

We start by shooting out a “solution” with \overline{x}_{1} = 0, \overline{x}_{2} = 0. This “solution” is computed to satisfy as many equations of the system as possible, while we move down the rows. From the second row, we get \overline{x}_{3} = a_{2}\overline{x}_{2} + b_{2}\overline{x}_{1} + c_{2}. We continue to solve for \overline{x}_{4}, ⋯ \overline{x}_{n}, recursively, by going down the rows. As we come to the (n − 1)st row, the “solution” \overline{x} = \{\overline{x}_{1}, \overline{x}_{2}, ⋯, \overline{x}_{n}\} has already been completely determined. As we plug the values of \overline{x}_{n−1} and \overline{x}_{n} into the last and then the first rows, chances are the two equations are not satisfied. The errors

ϵ_{1} = c_{n} + a_{n}\overline{x}_{n} + b_{n}\overline{x}_{n−1}  (3.2)

and

ϵ_{2} = c_{1} + b_{1}\overline{x}_{n}  (3.3)

indicate by how much we have missed the target. In other words, \overline{x} is not a correct solution to (3.1), but rather to

K\overline{x} = c + (ϵ_{2}, 0, ⋯, ϵ_{1})′.  (3.4)

Next we use the linearity of the operator K to seek a correction to the “solution” by considering the homogeneous system associated with (3.1):

Ky = 0.  (3.5)

We do a similar shooting for (3.5), with initial values \overline{y}_{1} = 1 and \overline{y}_{2} = 0. Again, we end up with a “solution” \overline{y} having errors ρ_{11}, ϱ_{21} ≠ 0 (otherwise, K would be singular):

L\overline{y} = (ρ_{21}, 0, ⋯, ρ_{11})′.  (3.6)

We do one more shooting with initial values \overline{z}_{1} = 0 and \overline{z}_{2} = 1 and obtain the corresponding errors ρ_{12} and ρ_{22}. The two error vectors form a matrix

R = ( \begin{matrix}ρ_{11} & ρ_{12} \\ ρ_{21} & ρ_{22} \end{matrix} ) .  (3.7)

We solve the 2 × 2 matrix equation

R ( \begin{matrix}α_{1} \\ α_{2} \end{matrix} ) = ( \begin{matrix}ϵ_{1} \\ ϵ_{2} \end{matrix} )  (3.8)

to find the correction needed to adjust our first shooting to hit the target. In fact, the linear combination

\overline{x} − α_{1}\overline{y} − α_{2}\overline{z}  (3.9)

is now a solution of (3.1).

It is easy to see how this method can, in theory at least, be extended to handle five-point stencils exemplified by the discrete Ginzburg-Landau operator. It can be verified that the discrete Ginzburg-Landau operator is a Hermitian operator on the complex vector space of dimension N × N (over N^{2} grid points). It is also not difficult to see that it is a negative definite operator. It is, therefore, nonsingular.

We first sweep (shoot) the nonhomogeneous system, starting from the right edge of Ω. As initial starting values, we need values of ψ to be given for apprpriate initial grid points. The most obvious choice is to let ψ = 0 for all points in the first two columns. Since each point in the second column is related to its four neighbors by the stencil and the only one that is not known at the moment is the right-hand neighbor, the latter can be easily determined by using the coefficients of the stencil and the nonhomogeneous part of the system of equations. The values of ψ on points in the third column are, therefore, determinable. We continue this sweeping action until the right boundary of Ω is reached, and wrap around to the first two columns to compute the errors. Note that the two columns give a total of 2N errors, to be combined to form a vector ϵ of length 2N.

Another choice of initial grid points is to take every other point in the first four columns in a checkerboard pattern. Every one of the remaining points in the second and third columns is surrounded by four initial points, and so the stencil can be used to determine ψ at that point.

To compute the rectifying matrix R, we go to the homogeneous stencil. For each grid point in the two initial columns, we do a sweep with initial value 1 at the the chosen point and 0 elsewhere. Each sweeping yields an appropriate column in R. As before, −R^{−1}ϵ gives the correct initial values that would give the exact solution to the system.

How much computation is involved in the algorithm? The initial sweeping of the nonhomogeneous system requires 4N multiplications and just as many additions. Each sweeping of the homogeneous system needs somewhat less, since on part of the columns in the left half of Ω, ψ will be zero, and nonhomogeneous terms are not present. The total is on the order of O(N^{2}). The main bulk of the computation is, therefore, in the inversion of R, which is on the order of O(N^{3}).

Parallelism can be exploited in many parts of the method. In the step for getting each new column in the sweeping, the computation for each point in the column can be done independently of the others. Then the 2N sweepings of the homogeneous system can be carried out in parallel.

A modification of the method allows the doubling of parallelism. One can sweep simultaneously in both directions on either side of the initial columns (provided that the coefficients of the stencil are nonzero, as is the case in the discrete Ginzburg-Landau operator). One starts with two columns at the center of the grid and sweeps both left and right; the error vector is computed as the sweeping fronts wrap around and meet.

One can also extend the well-known idea of multiple shooting to multiple sweeping. One starts with k sets of two initial columns spaced uniformly apart and sweeps in both directions. There are then 2kN unknowns to be solved.

The inverse of the five-point stencil is embedded in the inverse of R. In a numerical application in which the inverse of the operator is needed for more than a few times, it is more worthwhile to have R^{−1} computed explicitly and stored. This is the case in the solution of the Ginzburg-Landau model in which a Newton method is necessitated by the nonlinearity. In one approach, in which the two Ginzburg-Landau equations are to be solved simultaneously in each Newton cycle, the vector potential A, and hence the Ginzburg-Landau operator, changes from one step to another. In other words, the inverse R^{−1} has to be computed for each step. In practice, however, A does not change much from one step to another. If one uses the quasi-Newton method, one can use the same R^{−1} for several steps before computing a new one. In another approach, the Gauss-Seidel principle is used. Each of the two Ginzburg-Landau equations is solved in turn, for ψ and A, respectively, assuming the value of the other variable as found in the last iterative step. Then the same R^{−1} is used many times in each of the steps involving the solution for ψ.

If the coefficients of the stencil allow it, one can also start with two initial rows of grid points and sweep in the upward/downward direction.

In the numerical simulation of a time-dependent version of the Ginzburg-Landau equations, it is important to be able to find the largest eigenvalue of the discrete Ginzburg-Landau operator and the associated eigenvector. In a well-known algorithm for computing the smallest eigenvalue, the inverse of the operator is applied to a fixed vector a large number of times. Once R^{−1} is known, this algorithm can be applied to the Ginzburg-Landau operator without difficulty.

4 Partial Sweeping and Iterative Algorithms

We have tested the sweeping method using MATLAB and Fortran codes. In most cases, single-direction sweeping works satisfactorily for N up to about 12, and two-direction sweeping works for N up to about 18. Beyond this, the rectifying matrix R is so illconditioned that the inverse obtained is unreliable. As the columns are solved one by one, the test “solution” \overline{x} grows exponentially, giving the matrix R very large entries. The ratio between the largest and smallest eigenvalues of R can turn out to be extremely large. Details of the numerical experiments will be given in a forthcoming report.

As an example, we used the single-direction method and the first two columns as initial columns to solve the Ginzburg-Landau operator given by (see (2.6)–(2.12))

A(x, 0) = 0, A(x, L_{y}) = 2.25717254215814,  (4.1)
A is linear,  (4.2)
B(x, y) = 0,  (4.3)

and

n = 2, L_{y} = 3, N = 15.  (4.4)

The entry in R with the maximum absolute value is on the order of 7×10^{15}, and the ratio of the largest to the smallest eigenvalue of R is on the order of 10^{16}. Indeed, any kind of numerical computation on such an illconditioned matrix is highly unreliable. We found that even the associative law is not satisfied (numerically); MATLAB gives

max(abs((RO)R − R(OR)) ≈ 7 × 10^{16},  (4.5)

where O is the square matrix with all ones.

The sweeping method works better for an important class of five-point stencils, the Dirichlet stencils, which arise in partial sweeping algorithms discussed below. These are stencils for which the first column of L and the last column of R are zeros. Numerical experiments can handle well such stencils with N up to 40, by using two-direction, two-stage sweeping algorithms. The use of multi-stage algorithms will definitely permit larger N.

Although straightforward sweeping proves to be impractical for large systems, many ideas can be borrowed from other methods to obtain usable modifications. These include the technique of divide and conquer (in the form of domain decomposition) and alternating directions. Most of the procedures suggested below have been tried in sample cases. A thorough investigation will constitute an extensive future project.

In an iterative scheme, one computes a sequence of approximate solutions x_{n} → x (converging to the exact solution x ):

Lψ_{n} = c_{n} → c,  (4.6)

where L stands for the discrete Ginzburg-Landau operator or some suitable five-point stencil. In each step, instead of computing ψ_{n} directly, we look for a correction δ to be added to ψ_{n−1}, by solving

Lδ = c − Lψ_{n−1} = β_{n}  (4.7)

as accurately as possible. This is usually achieved by using the best approximate inverse for L found so far. There is the possibility that the error introduced in solving (4.7) may actually increase the error term β_{n+1} for the next step and that the algorithm will fail to converge. Suppose this is not the case; then the process is terminated when β_{n} becomes sufficiently small.

The method of domain decomposition implies a block decomposition of the matrices ψ_{n}, β_{n}, and others. The blocks need not be square, and they do not necessarily have equal dimensions. For instance, an N × N matrix can be decomposed into horizontal strips of dimension m×N or square blocks of dimension m × m, or a mixture of the two types. The idea is then to solve, for each block, only the subset of equations given by the stencil and only for those unknowns that fall within the subdomain. Referring to equation (4.7), we are, in fact, assuming the components of δ to be 0 outside the subdomain and attempt to change the values of those components inside the subdomain so as to satisfy as many equations as possible (all those that fall within the domain). In doing so, the stencil equations at all the grid points touching the subdomain may be violated. We describe this process as sweeping the errors represented by β_{n} in the subdomain to all the adjoining grid points. The hope, of course, is that the process of sweeping does not merely move the errors around, but reduces them through the adjustment of δ inside the subdomain.

Fortunately, for operators derivable from an optimization problem, such as the discrete Ginzburg-Landau operator, this hope is substantiated. Solving the stencil in a subdomain is equivalent to minimizing the energy functional under the constraint that the values of the unknown function are fixed at points outside the subdomain. Hence, after sweeping each subdomain, the functional is strictly decreased, and the total error is, in some appropriate sense, strictly decreased. It then follows that the error can only converge to zero. The idea of the method is definitely not new. It can be traced back, for example, to the method of balayage of Poincaré (see [6]) in proving the existence of a solution to the Laplace equation. A theoretical proof of convergence is, of course, not a guarantee for convergence in a numerical implementation; as usual, rounding errors can cause trouble. A typical safeguard is to have the functional computed after each iterative step and to make sure that it has really decreased.

Let us also give a heuristic justification for the method of partial sweeping. The fact that R is very illconditioned means that a slight change in the initial guesses, assigned on the initial columns, can produce, after sweeping, huge errors in the last column. To a lesser extent, a slight change in the values on the two columns adjoining the left or right boundaries of a subdomain corresponds to large errors at columns far away from the boundary, namely, those near the interior. Conversely, if one wishes to eliminate errors appearing in the interior of a subdomain, one does not expect much adjustment required to those points outside the subdomain. Hence, one may as well assume that they have fixed values when trying to reduce errors that occur in the interior of the subdomain. The argument, of course, fails for errors that are still within the subdomain but near the boundary. We shall describe below a possible way to deal with this situation.

In general, the convergence of the sequence of approximate solutions depends on the spectral radius of the partial operators that are being inverted and the relative errors incurred by the splitting of the operator.

The sweeping algorithm described in Section 3 has to be modified for partial sweeping. Let us first consider a subdomain bounded on all sides by grid points with fixed ψ values. The choice of sweeping direction is a matter of convenience, since all choices lead to identical results. Suppose we decide to sweep from the left. We now need only one initial column. A good choice is the second one next to the left boundary, since we can make use of the known values of ψ on the column just left of the boundary. Let m be the length of the columns. The last column at the right-hand boundary is used to give the error vector and the columns in the rectifying matrix R. The sweeping process, therefore, involves m unknowns and the inversion of the m × m matrix R. When compared with complete sweeping, we have fewer columns to sweep through (and so less room for the computed numbers to grow) and significantly fewer unknowns to solve for. Of course, we pay the price by having to iterate.

We can also do a two-direction sweeping, namely, inward from the left and the right boundaries of the subdomain. This will require solving for 2N unknowns, to give the initial values to be used at the first and the last columns. This can be done in parallel, though. The most important advantage of a two-direction sweeping is its increased stability.

We have also devised a two-stage sweeping in which one half of the domain is swept first and then the second half is swept without residual error in a way dependent on the outcome of the first sweeping. This is not an interative procedure. The idea can be easily extended to get a multi-stage sweeping.

The case of a subdomain that spans the entire length or width of ψ, such as a horizontal strip of size m × N, is slightly different. We have a choice of sweeping horizontally (and need two initial columns), leading to the inversion of a possibly illconditioned 2m × 2m matrix, or sweeping vertically, leading to the inversion of an N × N matrix. When N ≤ 2m, the obvious choice is the latter. When N is much larger than m, one may be tempted to use horizontal sweeping to avoid the excessive computation required in dealing with a large N × N matrix. However, one can see with a simple example that the N ×N matrix R is a banded matrix, allowing some known iterative methods to be employed. One can, in fact, apply another domain decomposition technique to the strip itself.

One very attractive feature of domain decomposition is that the subdomains can be swept independently of each other, and hence, in parallel. Thus, one iterative cycle consists of the simultaneous sweeping of all the subdomains in one decomposition.

We do not have to use the same decomposition in every cycle. In fact, this is not even advisable. After each cycle, the errors have been swept from the interior of each subdomain towards the boundary. Sometimes, the resulting errors may appear to be large, especially when compared to an approximate solution for which the errors are more or less evenly distributed among all the grid points. As we have argued above, errors swept from deep inside a domain have a better chance of being reduced during the process than errors that are present near the boundary. We can make use of this observation as we progress. The strategy we suggest is that in the next cycle, these boundary points (or as many of them as possible) be made the interior of the subdomains of a new decomposition.

As a concrete example, we can use two alternating decompositions, the first a familiar checkerboard pattern of m × m square blocks (where m is a factor of N), and the second obtained by shifting the square blocks by m/2 rows and m/2 columns. The half-size rectangles appearing on the top, as well as those on the left, are to be joined with the corresponding half on the bottom and the right respectively. Likewise, the four small squares at the four corners are to be united as one subdomain. As another example, the first decomposition is made of horizontal strips of size m × N and the second decomposition obtained by a shifting of m/2 rows. An alternative is to choose the second decomposition using vertical strips of size N ×m. This is an extension of the familiar method of alternating directions.

The idea can be extended: We can use multiple domain decomposition within a single computational cycle. This strategy is proposed to overcome the fact that partial sweeping is less effective for errors occurring near the boundaries of subdomains. Since the operator L is linear, one can break (4.7) up into a linear combination of similar systems, by breaking β_{n} up into several terms, each having nonzero entries only at the interiors of the subdomains of a suitable domain decomposition. All of these are then solved in parallel (and each with a parallel partial sweeping).

Since errors in the interior of a domain have a diminished effect on the required adjustments at the boundary grid points, we use yet another modification of the sweeping algorithm. Instead of computing R^{−1} and using it to determine the initial column, we compute matrices to give two or more interior columns. The remaining columns can be filled in by backward sweeping. We do not rely completely on backward sweeping because of its instability.

When the inversion of the original stencil is used in conjunction with Newton’s method, an exact inversion of the stencil is probably not necessary. Some errors are likely to arise in each Newton step anyway, so a moderately accurate inversion is all that is needed for practical purposes. This means that a few cycles of partial sweepings are usually sufficient in each Newton iteration.

Finally, we mention the idea of multigrid. Although we have not implemented this directly into our sweeping algorithm, it has been used in connection with solving the Ginzburg-Landau equations. The idea is to start with a coarser grid discretization to obtain a first approximation to the solution of the Ginzburg-Landau equation. A finer grid is then used, and the new Newton procedure starts with an initial guess given by interpolating the first approximate solution. Having a good initial guess greatly speeds up the convergence.

Acknowledgments I thank my colleagues Paul Plassmann and Steve Wright for many useful discussions and for sharing their recent work, especially in making available their numerical codes. I also thank Shirin R. Bar-Sela for writing some of the MATLAB programs and performing many of the numerical experiments. Miss Bar-Sela is a participant in the spring 1992 Student Research Participation Program at Argonne National Laboratory and is currently an undergraduate at the University of Houston – Downtown.

References

[1] Abrikosov, A., On the magnetic properties of superconductors of the second type, Zh. Eksperim. i Teor. Fiz., 32 (1957), 1442–1452. [English

translation: Soviet Phys. – JETP, 5 (1957), 1174–1182.]

[2] Bardeen, J., Cooper, L., and Schrieffer, J., Theory of superconductivity, Phys. Rev., 108 (1957), 1175–1204.

[3] Ascher, U. M., Mattheij, R. M. M., and Russell, R. D., Numerical Solutions of Boundary Value Problems for Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, New Jersey, 1991.

[4] Adler, S. L., and Piran, T., Relaxation methods for gauge field equilibrium equations, Reviews of Modern Phys., 56 (1984), 1–40.

[5] Chapman, S. J., Howison, S. D., and Ockendon, J. R., Macroscopic models for superconductivity, SIAM Review, to appear.

[6] Courant, R., and Hilbert, D. Methods of Mathematical Physics, Vol. II, Interscience Publishers, New York, 1962.

[7] DeGennes, P., Superconductivity in Metals and Alloys, Benjamin, New York, 1966.

[8] Doria, M., Gubernatis, J., and Ranier, D., Soving the Ginzburg-Landau equations by simulated annealing, Phys. Rev. B, 41 (1990), 6335–6340.

[9] Du, Q., Gunzburger, M., and Peterson, J., Analysis and approximation of a periodic Ginzburg-Landau model for type-II superconductors, SIAM Review, to appear.

[10] Du, Q., Gunzburger, M., and Peterson, J., Modeling and analysis of a periodic Ginzburg-Landau model for type-II superconductors, preprint, 1992.

[11] Frahm, H., Ullah, S. and Dorsey, A. T., Flux dynamics and the growth of the superconducting phase, Phys. Rev. Letters, 66 (1991), 3067–3070.

[12] Garner, J., Spanbauer, M., Benedek, R., Strandburg, K., Wright, S., and Plassman, P., Critical fields of Josephson-coupled superconducting multilayers, Mathematics and Computer Science Division Preprint MCS-P281-1292, Argonne National Laboratory, Argonne, Ill., December 1991.

[13] Ginzburg, V., and Landau, L., On the theory of superconductivity, Zh. Eksperim. i Teor. Fiz., 20 (1950), 1064–1082. [English translation: Men of Physics: L. D. Landau, ed. D. ter Harr, Pergamon, Oxford (1965), 138–167.]

[14] Gor’kov, L., Microscopic derivation of the Ginzburg-Landau equations in the theory of superconductivity, Zh. Eksperim. i Teor. Fiz., 36 (1959), 1918–1923. [English translation: Soviet Phys. – JETP, 9 (1959), 1364– 1367.]

[15] Kaper, H. G., Mathematical models of superconductivity, Argonne Technical Memorandum ANL/MCS-TM-146, March, 1991.

[16] Kato, R., Enomoto, Y., and Maekawa, S., Computer simulations of dynamics of flux lines in type II superconductors, preprint, 1991.

[17] Kuper, C., An Introduction of the Theory of Superconductivity, Clarendon, Oxford, 1968.

[18] Tinkham, M., Introduction to Superconductivity, McGraw-Hill, New York, 1975.

[19] Wang, Z. D., and Hu, C. R., A numerical relaxation approach for solving the general Ginzburg-Landau equations for type-II superconductors, preprint, 1991.