Multi-Constrained Variational Problems of Nonlinear Eigenvalue

Type : New Formulations and Algorithms

Alexander Eydeland, Joel Spruck and Bruce Turkington University of Massachusetts at Amherst

1 Introduction

The present work was motivated by a reexamination of some fundamental constrained variational problems that arise in equilibrium theory in fluid dynamics and Magnetohydrodynamics. These problems lead formally to nonlinear eigenvalue problems as variational equations. Indeed, it is traditional that many fundamental problems of mathematical physics can be formulated as nonlinear eigenvalue problems of the form

−Δu = Λ(u), u ∈ H_{0}^{1}(Ω)

The profile function Λ(u) is usually expressed in the form Λ(u) = λf(u) where f(u) is a known function and λ is a Lagrange multiplier.

However, these models do not always represent well the original physical problem. For in the original variational formulation, the constraints of the equilibrium formulation are derived from the conserved quantities associated to the evolution of the dynamical equations and thus the profile function is completely unknown and must be determined implicitly from the constraints.In MHD theory,such an objection was raised by Grad who introduced so called generalized or queer differential equation involving differentiations of the distribution function of u. These equations are purely formal unless u happens to be smooth with simple level sets. Nevertheless, Grad and his colleagues introduced numerical techniques to solve these equations.

As a prototypical problem of the type that that arises in many physical problems we pose the variational problem (P_{∞})

(1.1) E(u) = frac12 ∫_{Ω}(|∇u|^{2} dx → min over u ∈ M_{∞}(ū).

where for a given ū ∈ H_{0}^{1}(Ω), M_{∞}(ū) is the class of functions which are equimeasurable with ū. For simplicity we will assume thatū is non-negative.

The existence of a weak solution u ∈ H_{0}^{1}(Ω) of (P_{∞}) is straightforward . On the other hand, the construction of an appropriate variational equation (that is a Lagrange multiplier rule to determine the profile function Λ(u)) satisfied by such a solution is unclear. Also unclear is the question of the optimal regularity of the solution. The main source of difficulty in investigating these questions lies in the nature of the family of constraints. Our approach to this difficulty is to formulate a clear and conceptually simple classical variational problem (P_{n}) which replaces the family of constraints in (P_{∞}) by a finite family of simpler constraints, and which approximates (P_{∞}) when n is large. Besides being more tractable analytically, the problem (P_{n}) is also amenable to a numerical method of solution, the development of which is a principal goal in this paper. One of the key points of our approach is the construction of a finite element or constraint space in which the basis functions have very special properties.

We introduce a partition {σ_{i}}_{i=0}^{n} of the interval 0 ≤ σ < +∞ into n subintervals 0 = σ_{0} < σ_{1} < … < σ_{n−1} < σ_{n} = + ∞. For a partition of this kind, we introduce the constraint functionals

F_{i}(u): = ∫_{Ω} f_{i}(u)dx (i = 1, …, n) where (1.2)

f_{i}(s): = ∫_{σ_{i−1}}^{σ_{i}}(s − σ)_{+} dσ = l/2 (s − σ_{i−1})_{+}^{2} − l/2 (s − σ_{i})_{+}^{2} (1.3)

for s ∈ (we write s_{+} = max{s, 0} and s_{+}^{2} = (s_{+})^{2} .) Then we let

M_{n}(ū) = {u ∈ H_{0}^{1}(Ω): F_{i}(u) = F_{i}(ū) for all i = 1, …, n}. (1.4)

The problem (P_{n}) is the following multiconstrained minimization problem which may be viewed as a natural discretization of (P_{∞})

E(u) → min over u ∈ M_{n}(ū). (1.5)

Equivalently, (P_{n}) may be expressed as

(1.6) E(u) → min subject to F_{i}(u) = γ_{i} (i = 1, …, n),

where γ_{i}: = F_{i}(ū) = ∫_{σ_{i−1}}^{σ_{i}} β(σ) dσ may be viewed as given data. The existence of a solution of (P_{n}) and the form of the variational equation that it satisfies are given in the following result.

If ū is given so that γ_{i} > 0 for every i, then there exists a non-negative minimizer u ∈ C^{2,β} for (P_{n}) satisfying

−Δu = ∑_{i=1}^{n} λ_{i}fi^{′}(u) for some λ_{i} ∈ .

One point that should be emphasized is that our finite constraint problem yields exact solution of the relevant physical governing equations(for example, the ideal MHD equations) and is an approximation only in the sense that the full constraint space is not prescribed. Our approach is very different to that of Laurence and Stredulinisky who introduce an approximation to the Grad G.D.E. model that doesn’t enjoy this property.

As we have mentioned earlier, one of the principal goals of this paper is to develop an efficient and justifiable numerical method to compute the solutions of problem (P_{n}). To this end, we introduce an iterative procedure which is designed to converge to the solutions of P_{n}. Our algorithm A_{n} may be succinctly defined as follows:

Given u^{0} ∈ M_{n}(ū), u^{k+1} is the unique solution of the variational inequality

(1.7) E_{τ}(u) = E(u) + \frac{τ}{2}⟨u, u⟩ → min subject to

F_{i}(u^{k}) + ⟨F′_{i}(u^{k}), u − u^{k}⟩ ≥ γ_{i}, i = 1, …, n (1.8)

which is defined in terms of a positive parameter τ and with the admissible functions u lying in the constraint space M_{n}(ū) defined by the n linearized inequality constraints.

For a discussion of the main ideas underlying our construction of algorithm A_{n} the reader is referred to Appendix2. We do want to emphasize here that the algorithm is designed to take maximal advantage of the convexity of the basis functions f_{i}(s) even though the variational problem P_{n} is highly non-convex.

Perhaps the main result of the paper is the following

Given any initialization u^{0} ∈ M_{n}(ū), there exists a constant C = C(Ω, E(u^{0}), γ_{1}^{−1} … γ_{n}^{−1}, (Δσ)^{−1}) so that if τ > C, the sequence u^{k} generated by the algorithm A_{n} converges strongly in H_{0}^{1}(Ω) to the set S of critical points of the variational problem P_{n}, in the sense that dist(u^{k}, S) → 0 as k → ∞ where dist(u, S) = inf {||u − v||_{H_{0}^{1}(Ω)}: v ∈ S})

We now discuss the organization of the paper. In Sections 2-4, we restrict our discussion to a prototypical case of the class of variational problems that hold our attention throughout the paper. In Section 2 we formulate the prototype finite constraint problem (P_{n}), and establish its relation to the problem(P_{∞}. . Sections 2 contains the proof of Theorem 1.1 and illustrates in a very clear and concrete way, some of the fundamental ideas of classical optimization and duality theory. In particular, Lemma 2.3 and its relation to the linear independence of the constraints is of independent interest. In Section 3 we describe the algorithm (A_{n})for solving (P_{n}), which is a globally convergent iterative procedure. The basic convergence properties of this algorithm are derived and a more concrete description of the algorithm is given in terms of the so called dual variational problem. Section 4 contains the proof of the convergence of the algorithm to a solution of our original problem. The key step in the proof of the convergence is an apriori estimate (Lemma 4.3) which is of independent interest. In Section 5, we briefly indicate how to extend the development of Sections 2-4 to include a free boundary interface, to apply to distributions which change sign, and finally to apply to a more general class of variational problems. Section 6(Appendix 2)contains a brief description of the fundamental variational problems in MHD equilibrium theory which partially motivated our prototype problems. In future work we hope to elaborate on the method and ideas of the present paper to study the more general problems of Appendix 1. In Section 7(Appendix 1) we collect some remarks concerning (A_{n}) which help to clarify the construction of this algorithm.

The following conventional notations will be used in the sequel. Let Ω ⊆ ^{N}(N = 1, 2, …) be a bounded open set whose boundary, ∂Ω, is smooth enough (say C^{1+β}). We denote as usual Lebesgue measure on Ω by dx = dx_{1} … dx_{N}, x = (x_{1}, …, x_{N}) ∈ Ω, and for any measurable subset A ⊆ Ω we write |A| = ∫_{A} dx for its measure. The inner product and norm on L_{2}(Ω) are denoted by

⟨u, v⟩ = ∫_{Ω} u(x)v(x)dx , ||u|| = ⟨u, u⟩^{\frac{1}{2}}

The Sobolev space H^{m}(Ω)(m = 1, 2, …) is identified as the space of those functions in L^{2}(Ω) whose (weak) partial derivatives of order ≤ m are in L^{2}(Ω). The space H_{0}^{1}(Ω) consists of those functions which vanish in the weak sense on ∂Ω, and H^{−1}(Ω) is its dual with respect to the pairing ⟨., .⟩. The class of continuous functions on Ω̄ = Ω ∪ ∂Ω is denoted by C^{0}(Ω̄), and C^{m}(Ω̄) ⊂ C^{0}(Ω̄)(m = 1, 2, …) denotes the class of m times continuously differentiable functions. Also, C^{m,a}(Ω̄)(0 < α ≤ 1) is the subclass of C^{m}(Ω̄) consisting of those functions whose derivatives of order = m are Hölder continuous with exponent α; thus, C^{0,1}(Ω̄) denotes the class of Lipschitz continuous functions on Ω̄. All of these notations are adopted from the standard reference [?], where the reader is referred for further details, if necessary.

2 Formulation of problems (P_{n}) and (P_{∞}).

Let Ω be a domain in ^{n}. We begin by posing the problem (P_{∞}), which we adopt as a model for the type of variational problems presented in Appendix 1 (where the reader is referred for an indication of the source of these problems in mathematical physics). For a given function ū ∈ H_{0}^{1}(Ω), we consider the class of functions which are equimeasurable with ū; namely, we let

M_{∞}(ū) = {u ∈ H_{0}^{1}(Ω): | {u > σ} |=| {ū > σ} | for all σ ∈ }. (2.1)

The class M_{∞}(ū) consists of all functions H_{0}^{1}(Ω) which are rearrangements of ū (with respect to Lebesgue measure dx on Ω). This characterization can also be phrased in terms of the standard rearrangement mapping [?] which is u → u_{∗} where

u_{∗}(α): = inf {σ ∈ : | {u > σ} |> α} (0 ≤ α ≤ |Ω|). (2.2)

Then we see that

M_{∞}(ū) = {u ∈ H_{0}^{1}(Ω): u_{∗}(α) = ū_{∗}(α) for 0 ≤ α ≤ |Ω|}. (2.3)

For the purposes of the present study, however, we prefer to define the class M_{∞}(ū) in terms of a certain infinity family of integral constraints.

Lemma 2.1 The following statements are equivalent:

(a) u ∈ M_{∞}(ū)
(2.4)  (b) ∫_{Ω} ϕ(u)dx = ∫_{Ω} ϕ (ū) dx for all 0 ≤ ϕ ∈ C_{0}()
(c) ∫_{Ω}(u − σ)_{+} dx = ∫_{Ω} ū − σ)_{+} dx for all σ ∈ .

Proof. That (a) implies (b) is standard; (b) is equivalent to (a) when functions of the form ϕ(s) = {_{1 s>σ}^{0 s≤σ} are used, and so (b) holds for any function ϕ(s) which is an increasing limit of step functions (in fact, for any nonnegative Borel measurable function ϕ(s)). That (b) implies (c) is trivial. That (c) implies (a) follows from differentiation of (c) with respect to the variable σ. Precisely, the right derivative with respect to σ yields the desired result since

t^{−1}[(s − σ − t)_{+} − (s − σ)_{+}] ↑ {_{1 ≤σ}^{0 s≤σ} as τ ↓ 0};

then, by virtue of monotone convergence, it follows that

lim_{t→0^{+}} t^{−1}{∫_{Ω}(u − σ − t)_{+}dx − ∫_{Ω}(u − σ)_{+}dx} = |{u > σ}|

for every σ, and similarly with ū replacing u. This proves the claimed equivalence. □

The problem (P_{∞}) is a constrained minimization problem whose objective functional is E(u): = 1/2 ∫_{Ω} |∇u|^{2} dx, (2.5)

and whose family of constraints is defined by (2.4). For a given ū ∈ H_{0}^{1}(Ω) with ū ≥ 0 a.e. in Ω, we pose the variational problem (P_{∞})

E(u) → min over u ∈ M_{∞}(ū). (2.6)

We assume that ū is nonnegative for the sake of simplicity in the exposition here, deferring the general case until §4. Equivalently, (P_{∞}) may be written explicitly in terms of its constraints as (P_{∞})

(2.7) E(u) → min subject to ∫_{Ω}(u − σ)_{+} dx = β(σ), 0 ≤ σ < + ∞,

where β(σ): = ∫_{Ω}(ū − σ)_{+} dx may be viewed as given data.

The existence of a solution u ∈ H_{0}^{1}(Ω) of (P_{∞}) is straightforward (see Theorem 2. ). On the other hand, the construction of an appropriate variational equation (such as a Lagrange multiplier rule) satisfied by such a solution is not routine, and indeed remains an unanswered question. Equally unclear is the question of the (optional) regularity of the solution. The main source of difficulty in investigating these questions is, of course, the nature of the family of constraints. With this in mind, we therefore proceed to formulate the problem (P_{n}) which replaces the family of constraints in (P_{∞}) by a finite family of simpler constraints, and which approximates (P_{∞}) when n is large. Besides being more tractable analytically, the problem (P_{n}) is also amenable to a numerical method of solution, the development of which is a principal goal in this paper. Thus, (P_{n}) may be viewed as the natural discretization of (P_{∞}).

We introduce a partition {σ_{i}}_{i=0}^{n} of the interval 0 ≤ σ < +∞ into n subintervals 0 = σ_{0} < σ_{1} < … < σ_{n−1} < σ_{n} = + ∞. For a partition of this kind, we introduce the constraint functionals

F_{i}(u): = ∫_{Ω} f_{i}(u)dx (i = 1, …, n) where (2.8)

f_{i}(s): = ∫_{σ_{i−1}}^{σ_{i}}(s − σ)_{+} dσ = l/2 (s − σ_{i−1})_{+}^{2} − l/2 (s − σ_{i})_{+}^{2} (2.9)

for s ∈ . (We write s_{+} = max{s, 0} and simply s_{+}^{2} = (s_{+})^{2} for the sake of brevity.) Then we let

M_{n}(ū) = {u ∈ H_{0}^{1}(Ω): F_{i}(u) = F_{i}(ū) for all i = 1, …, n}. (2.10)

It is easy to verify that M_{∞}(ū) ⊆ M_{n}(ū) for any n, and that if a sequence of partitions is taken so that max_{1≤i≤n−1} (σ_{i} − σ_{i−1}) → 0 and σ_{n−1} → + ∞ as n → + ∞, then M_{∞}(ū) = ∩_{n}M_{n}(ū), the class M_{n}(ū) being defined by the n-th partition in the sequence. (If ū ∈ C^{0}(Ω̄), say, then it is not necessary that σ_{n−1} → + ∞, but rather it suffices that σ_{n−2} ≤ sup_{Ω} ū < σ_{n−1}.)

The problem (P_{n}) is the following multiconstrained minimization problem

E(u) → min over u ∈ M_{n}(ū). (2.11)

Equivalently, (P_{n}) may be expressed as

(2.12) E(u) → min subject to F_{i}(u) = γ_{i} (i = 1, …, n),

where γ_{i}: = F_{i}(ū) = ∫_{σ_{i−1}}^{σ_{i}} β(σ) dσ may be viewed as given data. The existence of a solution of (P_{n}) and the form of the variational equation that it satisfies are given in the next theorem.

Theorem 2.2 If ū is given so that γ_{i} > 0 for every i, then there exists a minimizer u for (P_{n}) satisfying

(a) u ∈ H_{0}^{1}(Ω) ∩ C^{2,β}(Ω̄),
(2.13)  (b) u ≥ 0 in Ω,
(c) −Δu = ∑_{i=1}^{n} λ_{i}fi^{′}(u) for some λ_{i} ∈ .

Proof. We invoke a standard argument to establish the existence of u. Let u_{j} ∈ M_{n}(ū) be a minimizing sequence in the sense that

E(u_{j}) = inf ⟨E(ũ): ũ ∈ M_{n}(ū)⟩.

The sequence u_{j}, being bounded in H_{0}^{1}(Ω), has a subsequence which converges weakly in H_{0}^{1}(Ω) to a limit u. We write this subsequence as u_{j} (after reindexing) and we have

E(u) ≤ lim_{j→∞} E(u_{j}) = inf ⟨E(ũ): ũ ∈ M_{n}(n̄)⟩,

F_{i}(u) = lim_{j→∞} F_{i}(u_{j}) = γ_{i} for each i = 1, …, n,

using respectively the lower semicontinuity of E with respect to weak convergence in H^{1} and the continuity of each F_{i} with respect to strong convergence in L_{2}. Thus, the limit u ∈ H_{0}^{1}(Ω) solves (P_{n}).

Property (2.9b) follows from the observation that u_{j} can be replaced by |u_{j}| in the above argument without changing the conclusion. (In fact, it is not difficult to see that any solution of (2.9a,c) must be nonnegative.)

The variational equation (2.9c) satisfied by u is the standard Lagrange multiplier rule. As is shown in the reference text [?], it suffices that the functionals E, F_{i}: H_{0}^{1}(Ω) → be C^{1} and that the gradients F′_{1}(u), …, F′_{n}(u) be linearly independent. It is immediate from their definitions that the functionals involved are Fréchet differentiable with

⟨E′(u), v⟩ = ∫_{Ω} ∇u · ∇v dx, ⟨F′_{i}(u), v⟩ = ∫_{Ω} f′_{i}(u)v dx where

f′_{i}(s) = (s − σ_{i=1})_{+} − (s − σ_{i})_{+} (s ∈ ). (2.14)

Thus, interpreting these derivatives in the sense of distributions, we have

E′(u) = −Δu ∈ H^{−1}(Ω), F′_{i}(u) = f′_{i}(u) ∈ H_{0}^{1}(Ω). (2.15)

Moreover, it is now clear that E, F_{i} are C^{1} functionals. The crucial linear independence of the gradients f′_{i}(u) (i = 1, …, n) will be deduced from the following

Lemma 2.3 Let f(s) ∈ C^{1}[0, ∞) ∩ C_{p}^{2}(0, ∞) satisfy f(0) = f′(0) = 0, f′′(s) ≥ 0. Then for any 0 ≤ u ∈ H_{0}^{1}(Ω)

∫_{Ω} f(u) dx ≤ c|Ω|^{\frac{1}{2}} ||f′′||_{L^{2}(Ω)} ∫_{Ω} |∇u|^{2} dx if N = 2 (2.16)

∫_{Ω} f(u) dx ≤ c||f′′||_{L^{\frac{N}{2}} (Ω)} ∫_{Ω} |∇u|^{2} dx if N ≥ 2 (2.17)

Proof. Suppose n ≥ 2 > p ≥ \frac{2n}{n+4}. Then by convexity of f,

∫_{Ω} f(u) dx ≤ ∫_{Ω} uf′(u) dx ≤ ||u||_{L^{\frac{np}{np−n+p}} (Ω)} ||f′(u)||_{L^{\frac{np}{n−p}} (Ω)} (2.18)

and by the Sobolev inequality

\begin{matrix}||f′(u)||_{L^{\frac{np}{n−p}} (Ω)} ≤ c ||∇(f′(u)||_{L^{p}(Ω)} \\ ≤ c ||f′′(u)||_{L^{\frac{2p}{2−p}} (Ω)} ||∇u||_{L^{2}(Ω)} \end{matrix} (2.19)

If N = 2 we take p = 1 and obtain from (2.18) (2.19)

\begin{matrix}∫_{Ω} f(u) dx ≤ c||u||_{L^{2}(Ω)} ||f′′(u)||_{L^{2}(Ω)} ||∇u||_{L^{2}(Ω)} \\ ≤ c||f′′(u)||_{L^{2}(Ω)} ∫_{Ω} |∇u|^{2} dx \end{matrix}

proving (2.16). If N > 2, we choose p = \frac{2N}{N+4} and obtain from (2.18) (2.19),

\begin{matrix}∫_{Ω} f(u) dx ≤ c||u||_{L^{\frac{2N}{N−2}} (Ω)} ||f′′(u)||_{L^{\frac{N}{2}} (Ω)} ||∇u||_{L^{2}(Ω)} \\ ≤ c||f′′(u)||_{L^{\frac{N}{2}} (Ω)} ∫_{Ω} |∇u|^{2} dx \end{matrix}

□

Indeed, applying Lemma (2.3) to f_{i}(u) establishes that each of the sets {σ_{i−1} < u < σ_{i}} has strictly positive measure (with a lower bound depending only on upper bounds for γ_{1}^{−1}, …, γ_{n}^{−1}, E(u) and |Ω|).

Thus, if for some c_{i} ∈ the identity 0 = ∑_{i=1}^{n} c_{i}f′_{i}(u) holds a.e. in Ω , then by restricting their identity to the sets {σ_{i−1} < u < σ_{i}} we find successively that c_{1} = 0, c_{2} = 0, …, c_{n} = 0. Hence, the required linear independence of the gradients of the constraint functionals is proved, and so (2.9c) follows for some Lagrange multipliers λ_{i} (i = 1, …, n) uniquely determined by u.

Finally, the smoothness of u stated in (2.9a) can be derived from the standard regularity theory for elliptic partial differential equations. We introduce the Green function g(x, x′), x, x′ ∈ Ω, defined by −Δ_{x}g(x, x′) = γ(x − x′) for x ∈ Ω, g(x, x′) = 0 for x ∈ ∂Ω, then we have

u(x) = ∫_{Ω} g(x, x′)Λ(u(x′)) dx with Λ(s): = ∑_{i=1}^{n} λ_{i}f′_{i}(s). (2.20)

Since ||Λ(u)|| ≤ (max_{1≤i≤u} |λ_{i}|) ||u|| we obtain Λ(u) ∈ L^{2}(Ω), and hence u ∈ C^{0}(Ω̄), using elementary properties of g. Now differentiating (2.12) in x we find that u ∈ C^{1}(Ω̄). Thus, Λ(u) ∈ C^{0,1}(Ω̄), and this is optional since the Lipshitz function Λ(s) is piecewise linear. The statement (2.9a) that u ∈ C^{2,β}(Ω̄), is therefore a consequence of the global Schauder estimates applied to the equation (2.9c). (It is assumed here that ∂Ω is of class C^{2,β}, at least.) This completes the proof of (2.9). □

The function Λ(s) occurring in (2.12) will play a key role in our subsequent analysis, and will be referred to as the profile function associated with a solution u. For later reference we record its definition:

Λ(s) = \{ \begin{matrix}0 & if s < 0 \\ ∑_{i=1}^{n} λ_{i}f′_{i}(s) & if s ≥ 0, \end{matrix} (2.21)

where λ_{1} …, λ_{n} are real constants (the multipliers in Theorem 1.2). Such a profile function is piecewise linear with slope Λ′(s) = λ_{i} on each interval σ_{i=1} < s < σ_{i}. (We recall that the partition of [0, + ∞) into intervals [σ_{i−1}, σ_{i}), i = 1, …, n, is indexed so that 0 = σ_{0} < σ_{1} < … < σ_{n−1} < σ_{n} = + ∞.) Clearly, Λ(s) is nondecreasing on [0, + ∞) if and only if λ_{i} ≥ 0 for every i = 1, …, n. The monotone functions f′_{i} (s) (i = 1, …, n) may be viewed as forming a basis for the space of profile functions - namely, those functions Λ(s) which satisfy (i) Λ ∈ C^{0}(), (ii) Λ(s) = 0 for s ≤ 0, (iii) Λ(s) is linear (actually affine) on each interval σ_{i=1} ≤ s ≤ σ_{i}. This basis, each member of which is monotone, is related to the the standard basis consisting of finite-elements ϕ_{i}(s) for the partition (or grid) {σ_{i}}_{i=0}^{n} by the formulas

\{ \begin{matrix}ϕ(s) = f′_{i}(s)/Δσ_{i} − f′_{i+1}(s)/Δσ_{i+1} (i + 1, …, n − 2) \\ ϕ_{n−1}(s) = f′_{n−1}(s)/Δσ_{n−1}, ϕ_{n}(s) = f′_{n}(s), \end{matrix}

where Δσ_{i} = σ_{i} − σ_{i−1} denotes the increments. Then, at least when 1 ≤ i ≤ n−2, each piecewise linear ϕ_{i}(s) is supported on the interval σ_{i−1} ≤ s ≤ σ_{i+1} and is normalized by ϕ_{i}(σ_{i}) = 1. In terms of these finite-elements the profile function is represented as

Λ(s) = \{ \begin{matrix}0 & if s < 0 \\ ∑_{i=1}^{n−1} Λ_{i}ϕ_{i}(s) + λ_{n}ϕ_{n}(s) & if s ≥ 0, \end{matrix}

where Λ_{i} = Λ(σ_{i}); in other words, Λ(s) on 0 ≤ s ≤ σ_{n−1} is satisfied as the linear interpolant of its values Λ_{i} at the partition points (or grid nodes) σ_{i}.

The representation (2.13) in terms of the monotone basis functions f′_{i}(s) will be used throughout the sequel, however; its importance stems from the fact that the constraint functionals F_{i}, like the objective functional E, are convex for this choice of the functions f_{i}(s).

Furthermore, it is a consequence of the definition of f_{i}(s) in (2.7) that

∑_{i=1}^{n} f_{i}(s) = \frac{1}{2} s_{+}^{2} (s ∈ ). (2.22)

Therefore, any admissible function u ∈ M_{n}(ū) satisfies

\frac{1}{2} ∫_{Ω} u^{2}dx = ∑_{i=1}^{n} γ_{i} = \frac{1}{2} ∫_{Ω} ū^{2}dx. (2.23)

Indeed, if n = 1 (and so σ_{0} = 0 and σ_{1} = + ∞) then the class M_{1}(ū) is characterized precisely by the single constraint (2.15). In view of this fact, we may say that the problem (P_{n}) constitutes a generalization or extension of the classical Rayleigh principle [?] characterizing variationally the first eigenvalue-eigenfunction pair for the Laplacian operator −Δ on H_{0}^{1}(Ω). When n = 1, the problem (P_{1})

(2.24) ∫_{Ω} |∇u|^{2}dx → min subject to ∫_{Ω} u_{+}^{2} = 2γ_{1}

yields a solution pair (u^{(1)}, λ_{1}^{(1)}) ∈ H_{0}^{1}(Ω) × to the usual linear eigenvalue problem

−Δu^{(1)} = λ_{1}^{(1)})u^{(1)} in Ω, u^{(1)} = 0 on ∂Ω,

which, in addition, is nonnegative in Ω. When n ≥ 2, the problem (P_{n}) results in a solution pair (u^{(n)}, λ^{(n)}) ∈ H_{0}^{1}(Ω) ×^{n} of

(2.25) −Δu^{(n)} = Λ^{(n)}(u^{(n)}) in Ω, u^{(n)} = 0 on ∂Ω, u^{(n)} ≥ 0 in Ω,

a nonlinear eigenvalue problem for which the vector of eigenvalues λ^{(n)} = (λ_{1}^{(n)}, …, λ_{n}^{(n)}) determines the (piecewise linear) profile function Λ^{(n)}(s) according to (2.13). In going from the linear problem (P_{1}) to the general (and nonlinear) problem (P_{n}) the familiar L^{2}(Ω)-normalization ||u||^{2} = 2γ_{1} (=1, say) is replaced by the family of constraints F_{i}(u) = γ_{i}(i = 1, …, n), which in effect amounts to constraining n L^{2}(Ω)-expressions involving the lower functions (u − σ_{i})_{+} of u; namely, the constraints for (P_{n}) can be written as

||(u − σ_{i−1})_{+} ||^{2} − || (u − σ_{i})_{+} ||^{2} = 2γ_{i} (i = 1, …, n). (2.26)

The family of constraints (2.17) is therefore seen to be a strengthening of the single constraint (2.15) in the usual Rayleigh principle (P_{1}).

The sense in which the solutions u^{(n)} of (P_{n}) approximate solutions u^{(∞)} (say) of (P_{∞}) as n → + ∞ can now be determined. As mentioned above, we suppose that a sequence of partitions is taken with max_{1≤i≤n−1}(σ_{i−1} − σ_{i}) → 0 and σ_{n−1} → + ∞ as n → ∞ so that we are assured that M_{∞}(ū) = ∩_{n}M_{n}(ū). By virtue of the bound E(u^{(n)}) ≤ E(ū) < + ∞, every subsequence of {u^{(n)}}_{n=1}^{∞} has a further subsequence that converges weakly in H_{0}^{1} and strongly in L^{2}. If we let u^{(∞)} ∈ H_{0}^{1}(Ω) denote such a limit point of the solution sequence u^{(n)} for the multi-constrained problems (P_{n}) then we claim that u^{(∞}) is a solution of the (infinitely-constrained) problem (P_{∞}). Indeed, any admissible function ũ ∈ M_{∞}(ū) for (P_{∞}) is admissible for each problem (P_{n}), and hence, using the weak - H_{0}^{1} convergence, E(u^{(∞}) ≤ lim inf E(u^{(n)}) ≤ E(ũ), where n is understood to be tend to infinity along the subsequence. Also, using the strong -L^{2} convergence, it is straightforward to verify that u^{(∞)} ∈ M_{∞}(ū) because it is evident that for any σ ∈ (0, + ∞)

∫_{Ω}(u − σ)_{+} dx = lim_{n→∞} \frac{1}{Δσ_{i}} ∫_{Ω} f_{i}(u) dx,

where i = i(n, σ) is chosen so that σ ∈ [σ_{i−1}, σ_{i}) for each sufficiently large n. (We recall that Δσ_{i}: = σ_{i} − σ_{i−1}.) Thus the claim that u^{(∞)} solves (P_{∞}) follows. This result can be stated in the form

dist(u^{(n)}, S_{∞}) → 0 as n → + ∞, (2.27)

where S_{∞} = S_{∞}(ū) is defined to be the set of solutions of (P_{∞}) for a given function ū; the distance function dist(u, S) from a point u ∈ L^{2}(Ω) set S ⊆ L^{2}(Ω) is defined by

dist(u, S) = inf {||u − v||: v ∈ S}. (2.28)

That (2.18) is true is obvious from the above discussion. We remark that it is necessary to consider the set of solutions of (P_{∞}) since the uniqueness of its solutions has not been established and, in fact, may not hold in general. The statement (2.18) allows us to conclude, at least, that for large n the solution u^{(n)} of (P_{n}) approximates some solution of (P_{∞}) in the L^{2}(Ω)-norm. This conclusion justifies in part the discretization of the constraints for (P_{∞}) resulting in the finitely-many constraints for (P_{n}).

Before passing to the discussion of the algorithm for solving (P_{n}) we note that a variant of this problem can also be employed. We now assume that ū ∈ H_{0}^{1}(Ω) ∩ C^{0}(Ω̄) and let σ_{n} = sup_{Ω} ū (rather than σ_{n} = + ∞). With respect to a partition 0 = σ_{0} < σ_{1} < ⋯ < σ_{n−1} < σ_{n} < + ∞, we consider the problem (P′_{n})

E(u) → min over u ∈ M_{n}(ū), u ≤ σ_{n} a.e. in Ω.

This problem is of variational inequality type [?], being an obstacle problem with no additional constraints. It can be shown that the analogue of Theorem 1.2 holds in the sense that there is a solution u ∈ H_{0}^{1}(Ω) ∩ C^{1,1}(Ω̄) satisfying the equation

−Δu = ∑_{i=1}^{n} λ_{i}f′_{i}(u) in {u < σ_{n}} ⊆ Ω.

(If {u = σ_{n}} has positive measure, then u is not necessarily twice continuously differentiable, of course.) However, we shall restrict our further discussion to (P_{n}), leaving the parallel development for (P′_{n}) to the reader.

3 Construction of the algorithm for solving (P_{n}).

We now proceed to define an iterative procedure which is designed to converge to the solutions of (P_{n}). In this section we concentrate on describing the algorithm in its general form and establishing its basic convergence property. For a discussion of the main ideas underlying the construction of the algorithm the reader is referred to Appendix 2.

For a (fixed) positive constant τ, let E_{τ} be the following modified objective functional

E_{τ}(u) = E(u) + \frac{τ}{2} ||u||^{2} = \frac{1}{2} ∫_{Ω} [|∇u|^{2} + τu^{2}] dx. (3.1)

Also, for any v ∈ H_{0}^{1}(Ω) let L_{n}(v) be the following class of admissible functions

L_{n}(v) = {u ∈ H_{0}^{1}(Ω): F_{i}(v) + ⟨F′_{i}(v), u − v⟩ ≥ γ_{i} for all i}. (3.2)

We notice that in contrast to the class M_{n}(ū) defined by n nonlinear equality constraints, the class L_{n}(v) is defined by n linear (actually affine) inequality constraints. (L_{n}(v) depends upon ū through the constants γ_{i} = F_{i}(ū), but this dependence will not be emphasized in the notations.) In terms of these definitions we construct a sequence of approximations (u^{k}, λ^{k}) ∈ H_{0}^{1}(Ω) ×^{n} to the solutions of the variational problem (P_{n}) as follows: Algorithm (A_{n}) : given u^{0} ∈ M_{n}(ū), let (u^{k}, μ^{k}) ∈ H_{0}^{1}(Ω) × [0, + ∞)^{n} be defined iteratively by solving the sequence of convex optimization problems: u^{k+1} is the (unique) solution of

E_{τ}(u) → min over u ∈ L_{n}(u^{k}) (3.3)

and μ^{k+1} is the corresponding n-vector of nonnegative multipliers (uniquely determined by u^{k+1}); then let

λ^{k+1} = μ^{k+1} − τ. (3.4)

The algorithm (A_{n}) produces a well-defined sequence of approximations once an initialization u^{0} and a constant τ are chosen. If suffices to take u^{0} = ū. It will be shown in §4 that τ can be chosen sufficiently large to ensure the convergence of the iterative sequence (u^{k}, λ^{k}) to the set S of solution of the problem

\{ \begin{matrix}E′(u) = ∑_{j=1}^{n} λ_{j}F′_{j}(u) for some λ_{1}, …, λ_{n} ∈ \\ F_{i}(u) = γ_{i} (i = 1, …, n); \end{matrix} (3.5)

namely, the set of all critical points for the problem (P_{n}) which satisfy the given constraints (or, equivalently, belong to M_{n}(ū)). In the present development τ will be viewed as a (fixed) positive parameter.

The iterate (u^{k+1}, μ^{k+1}) in the algorithm (A_{n}) is characterized as the solution of the equations

\{ \begin{matrix}E′(u^{k+1}) + τu^{k+1} = ∑_{j=1}^{n} λ_{j}^{k+1}F′_{j}(u^{k}) \\ λ_{i}^{k+1} [F_{i}(u^{k}) + ⟨F′_{i}(u^{k}), u^{k}⟩ − γ_{i}] = 0 (i = 1, …, n) \end{matrix} (3.6)

Indeed, these equations are precisely the Kuhn-Tucker conditions associated with the convex optimization problem (2.22) which defines (A_{n}). The reader is referred to [?] for a proof that (2.22) and (2.24) are equivalent provided that there exists some ũ ∈ L_{n}(u^{k}) for which

F_{i}(u^{k}) + ⟨F′_{i}(u^{k}), ũ − u^{k}⟩ > γ_{i} for every i = 1, …, n.

This so-called Slater condition follows from the convexity of the functional F_{i}(u), which implies that

F_{i}(u^{k}) ≥ F_{i}(u^{k+1}) + ⟨F′_{i}(u^{k−1}), u^{k} − u^{k−1}⟩ ≥ γ_{i} (k = 1, 2, …). (3.7)

In addition, F_{i}(u^{k}) ≤ ⟨F′_{i}(u^{k}), u^{k}⟩ since F_{i}(0) = 0. thus, it is readily verified that the required condition is satisfied by ũ = (1 + ϵ)u^{k} with ϵ > 0. It should be emphasized that the equations (2.24) characterize the iterate (u^{k+1}, μ^{k+1}) ∈ L_{n}(u^{k}) × [0, + ∞)^{n}; and it is implicitly assumed in (2.24) that F_{i}(u^{k}) + ⟨F′_{i}(u^{k}), u^{k+1} − u^{k}⟩ ≥ γ_{i}, (3.8)

μ_{i}^{k+1} ≥ 0, (3.9)

for every i = 1, …, n. The vector μ^{k+1} = (μ_{1}^{k+1}, …, μ_{u}^{k+1}) is called the Kuhn-Tucker vector corresponding to the solution u^{k+1} of (2.22).

The above construction of the algorithm (A_{n}) depends fundamentally on the convexity of the optimization problem (2.22); namely, it relies on the strict convexity of the objective functional E_{τ} and the convexity of the class of admissible functions L_{n}(u^{k}), which is defined by n affine inequality constraints. The uniqueness of a minimizer u^{k+1} for (2.22) follows immediately from these attributes. In turn, the uniqueness of the Kuhn-Tucker vector μ^{k+1} follows from the smoothness of the objective and constraint functionals and the linear independence of the gradients F′_{i}(u^{k}). The latter property can be proved exactly as in the proof of Theorem 1.2, now using (2.25) with γ_{i} > 0.

We describe next a more explicit construction of the iterative procedure (A_{n}) which furnishes us with a concrete numerical implementation of the algorithm. Let G_{τ} denote the Green operator for the elliptic boundary value problem −Δw + τw = h in Ω, w = 0 on ∂Ω; (3.10)

that is, the solution is represented as w = G_{τ}h where G_{τ}: L^{2}(Ω) → H_{0}^{1}(Ω) ∩ H^{2}(Ω). For any v ∈ H_{0}^{1}(Ω) let

a_{ij}(v): = ⟨f′_{i}(v), G_{τ}f′_{j}(v)⟩ (3.11)

c_{i}(v): = γ_{i} − F_{i}(v) + ⟨F′_{i}(v), v⟩. (3.12)

for i, j = 1, …, n. (The dependence of these expressions on τ is left implicit, for the sake of simplicity in the notation.) We consider the quadratic form

Q(μ; v): = \frac{1}{2} ∑_{i,j=1}^{n} a_{ij}(v)μ_{i}μ_{j} − ∑_{i=1}^{n} c_{i}(v)μ_{i} (μ ∈ ^{n}) (3.13)

This defines a positive-definite quadratic form on^{n} whenever |{x ∈ Ω: σ_{i−1} < v(x) < σ_{i}}| > 0 for every i = 1, …, n. We check this fact by calculating

∑_{i,j=1}^{n} a_{ij}(v)μ_{i}μ_{j} = ⟨μ · f′(v), G_{τ}[μ · f′(v)]⟩

where we write μ · f′(v) = ∑_{i=1}^{n} μ_{i}f′_{i}(v). We then notice that the latter expression is zero if and only if μ · f′(v) = 0 a.e. in Ω. Arguing as in the proof of Theorem 1.2, we deduce that μ_{1} = 0, μ_{2} = 0, …, μ_{n} = 0 successively, as required.

We now claim that algorithm (A_{n}) is equivalent to the following explicit iterative procedure: for u^{0} ∈ M_{n}(ū), let (u^{k}, μ^{k}) ∈ H_{0}^{1}(Ω) × [0, + ∞)^{n} be defined iteratively by

u^{k+1} = ∑_{j=1}^{n} μ_{j}^{k+1} w_{j} (3.14)

where w_{j} = G_{τ} f′_{j}(u^{k}), (3.15)

\begin{matrix}μ_{j}^{k+1} ≥ 0 and Q(μ^{k+1}; u^{k}) \\ = min\{Q(μ; u^{k}): μ_{i} ≥ 0, i = 1, …, n\}. \end{matrix} (3.16)

In other words, (u^{k+1}, μ^{k+1}) is constructed from u^{k} in a three-step process:

(1) each w_{j}(j = 1, …, n) is found by solving (2.28) with h = f′_{j}(u^{k}); (2) μ^{k+1} is taken to be the unique solution of the quadratic programming problem (2.34); (3) u^{k+1} is assembled from μ_{j}^{k+1} and w_{j}(j = 1, …, n) according to (2.32). The verification of the claim that this procedure is equivalent to (A_{n}) as stated earlier is standard. Calculating the variational inequalities satisfied by the minimizer μ^{k+1} of (2.34) we have for each i = 1, …, n

∑_{j=1}^{n} a_{ij}(u^{k})μ_{j}^{k+1} − c_{i}(u^{k}) \{ \begin{matrix}≥ 0 & if μ_{i}^{k+1} = 0 \\ = 0 & if μ_{i}^{k+1} > 0. \end{matrix} (3.17)

Also, combining (2.32) and (2.33) we have

−Δu^{k+1} + τu^{k+1} = ∑_{j=1}^{n} μ_{j}^{k+1}f′_{j}(u^{k}) in Ω, u^{k+1} = 0 on ∂Ω. (3.18)

It is clear that (2.36) is a restatement of the first equation in (2.24); in turn, it is evident from the definitions of a_{ij} and c_{i} that (2.35) is equivalent to the second equation in (2.24) together with the inequalities (2.26) and (2.27) which supplement (2.24). Thus the claimed equivalence of the two forms of the algorithm is verified.

It is interesting to note that the optimization problems (2.34) and (2.22) are dual problems in the sense of convex analysis [?]. Indeed, we can directly verify that

−Q(μ; u^{k}) = min{E_{τ}(u) − ∑_{i=1}^{n} μ_{i}[⟨F′_{i}(u^{k}), u⟩ − c_{i}(u^{k})]: u ∈ H_{0}^{1}(Ω)}.

Then general theory informs us that the Kuhn-Tucker vector μ^{k+1} corresponding to the minimizer u^{k+1} for (2.22) is itself the maximizer for the dual problem:

−Q(μ; u^{k}) → max over μ_{i} ≥ 0 (i = 1, …, n).

The preceding construction of the explicit form of the iterative procedure (A_{n}) may therefore be viewed as a specific case of the general duality theory. We shall not need this level of generality in the sequel, however.

The algorithm (A_{n}) is devised to exploit as much as possible the important fact that the objective and constraint functionals for the variational problem (P_{n}) are convex. Even though (P_{n}) itself is not a convex optimization problem because its constraints are nonlinear equalities, its convexity attributes are enough to imply that the iterative sequence defined by (A_{n}) has very special monotonicity and convergence properties. These properties are the focus of our attention in the next two lemmas. Lemma 3.1 For every k we have

E_{τ}(u^{k+1} − u^{k}) ≤ E_{τ}(u^{k} − E_{τ}(u^{k+1}) (3.19)

E(u^{k+1}) ≤ E(u^{0}). (3.20)

Proof. To prove (2.37) we use the identity

E_{τ}(u^{k} − E_{τ}(u^{k+1}) = ⟨E′_{τ}(u^{k+1}), u^{k} − u^{k+1}⟩ + E_{τ}(u^{k+1} − u^{k}). (3.21)

By virtue of (2.25) we have u^{k} ∈ L_{n}(u^{k}). Therefore, (1 − t)u^{k+1} + tu^{k} ∈ L_{n}(u^{k}) for all 0 ≤ t ≤ 1, as this class of functions is convex. Consequently, since u^{k+1} solves (2.22) we get

\begin{matrix}E_{τ}(u^{k+1}) ≤ E_{t}(u^{k+1}) + t(u^{k} − u^{k+1})) = \\ E_{τ}(u^{k+1}) + t⟨E′_{τ}(u^{k+1}), u^{k} − u^{k+1}⟩ + 0(t^{2}) \end{matrix}

as t → 0+. Hence, ⟨E_{τ}(u^{k+1}), u^{k+1} − u^{k+1}⟩ ≥ 0, and so (2.37) follows from (2.39).

To prove (2.38) we note that

\frac{1}{2} ||u^{0}||^{2} = ∑_{i=1}^{n} γ_{i} ≤ ∑_{i=1}^{n} F_{i}(u^{k+1}) = \frac{1}{2} ||u^{k+1}||^{2}.

This yields the desired inequality since

E(u^{k+1}) = E_{τ}(u^{k+1}) − \frac{τ}{2} ||u^{k+1}||^{2} ≤ E_{τ}(u^{0}) − \frac{τ}{2} ||u^{0}||^{2} = E(u^{0}).

□

The above lemma establishes that for any initialization u^{0} ∈ M_{n}(ū) the sequence E_{τ}(u^{k}) is nonincreasing. This monotonicity property of the algorithm (A_{n}) leads directly to the following (partial) convergence property of the iterates u^{k}.

Lemma 3.2 Let u ∈ H_{0}^{1}(Ω) be any H_{0}^{1}-weak limit point of the sequence u^{k}. Then there exists a unique vector μ = (μ_{1}, …, μ_{n}) ∈ ^{n} such that

E′_{τ}(u) = ∑_{j=1}^{n} μ_{j} F′_{j}(u) (3.22)

μ_{i} [F_{i}(u) − γ_{i}] = 0 (i = 1, …, n) (3.23)

μ_{i} ≥ 0, F_{i}(u) ≥ γ_{i} (i + 1, … n). (3.24)

Proof. We recognize that (2.40), (2.41), (2.42) are precisely the Kuhn-Tucker conditions for the convex minimization problem

E_{τ}(ũ) → min over ũ ∈ L_{n}(u), (3.25)

and that μ is the Kuhn-Tucker vector corresponding to the minimizer u. Therefore, it suffices to show that u is indeed the solution of (2.43). Let u^{k_{p}} → u weakly in H_{0}^{1}(Ω), and hence strongly in L^{2}(Ω), as p → + ∞. Given an arbitrary ũ ∈ L_{n}(u) we consider the perturbation ũ + ϵu with ϵ > 0. Recalling the definition of L_{n}(u), we find that for each i = 1, …, n

F_{i}(u) + ⟨F′_{i}(u), (ũ + ϵu) − u⟩ ≥ γ_{i} + ϵ⟨F′_{i}(u), u⟩ ≥ (1 + ϵ)γ_{i}.

This implies that for sufficiently large p

F_{i}(u^{k_{p}}) + ⟨F′_{i}(u^{k_{p}}), (ũ + ϵu) − u^{k_{p}}⟩ ≥ (1 + ϵ/2)γ_{i},

by virtue of the continuity of the terms on the left-hand side of the latter inequality with respect to strong L^{2} convergence. Consequently, ũ + ϵu ∈ L_{n}(u^{k_{p}}) for sufficiently large p, and so we obtain

E_{τ}(u) ≤ E_{τ}(u^{k_{p}+1}) ≤ E_{τ}(ũ + ϵu),

where, in the first inequality, we invoke the monotonicity property (2.37) and the lower semi-continuity of E_{τ} with respect to weak H_{0}^{1} convergence. Since this inequality holds for arbitrarily small positive ϵ, we conclude that the limit point u solves (2.43), as required. □

The partial convergence result given in the above lemma provides the first step in the proof that the iterates (u^{k}, λ^{k}) defined by the algorithm (A_{n}) converge in an appropriate sense to the critical points of (P_{n}). Indeed, Lemmas 2.3 and 2.4 together permit us to conclude that any subsequence of the iterative sequence {u^{k}} ⊆ H_{0}^{1}(Ω) has a further subsequence which converges (weakly in H_{0}^{1} and strongly in L^{2}) to a solution u of (2.40), (2.41), (2.42). But (2.42) can be written equivalently as

E′(u) = ∑_{j=1}^{n} λ_{j} F′_{j}(u) where λ_{j}: = μ_{j} − τ, (3.26)

in view of the basic identity (2.14) which implies that

∑_{i=1}^{n} f′_{j}(s) = s_{+} (s ∈ ). (3.27)

(Here we also use the fact that u^{k} ≥ 0 for every k, and hence u ≥ 0, as follows immediately from the explicit form of (A_{n}) - namely, (2.32), (2.33), (2.34). Thus the addition of the term \frac{τ}{2}||u||^{2} to the objective functional E(u) has the effect of shifting the multipliers by τ, sending λ_{i} into μ_{i} = λ_{i}+τ. Accordingly, a solution u of (2.40), (2.41), (2.42) is a critical point of the variational problem (P_{n}) provided that μ_{i} > 0 for every i = 1, …, n, since then the constraints must be equalities F_{i}(u) = γ_{i} for every i = 1, …, n. The latter condition can be ensured by choosing τ large enough (τ > max {0, −λ_{1}, …, −λ_{n}}) depending on the multipliers γ_{i} corresponding to any solution u of (2.44). The a priori estimate on max_{1≤i≤n} |λ_{i}| needed to complete this argument is given in Lemma 3.2, and so it provides the second step in the proof of convergence of the algorithm (A_{n}). The complete convergence result is contained in Theorem 4.1.

4 Convergence Theorems

In the previous section we have studied the algorithm A_{n}:

Given u^{0} ∈ M_{n}(ū), u^{k+1} is the unique solution of the variational problem

\begin{matrix}E_{τ}(u) = \frac{1}{2} ∫_{Ω}(|∇u|^{2} + τu^{2}) → min \\ u ∈ H_{0}^{1}(Ω): F_{i}(u^{k}) + ⟨F′_{i}(u^{k}), u − u^{k}⟩ ≥ γ_{i}, i = 1, …, n \end{matrix} (4.1)

which is defined in terms of a positive parameter τ and with linearized inequality constraints.

In this section we will show that τ can be chosen apriori so that the limit points of the Algorithm A_{n} are solutions of our original problem P_{n}. More precisely, we have

Theorem 4.1 Given any initialization u^{0} ∈ M_{n}(ū), there exists C = C(Ω, E(u^{0}), γ_{1}^{−1} … γ_{n}^{−1}, (Δσ)^{−1}) so that if τ > C, the algorithm A_{n} given by (4.1) converges strongly in H_{0}^{1}(Ω) to the set S of critical points of our original problem P_{n}, in the sense

dist(u^{k}, S) → 0 as h → ∞
where dist(u, S) = inf {||u − v||_{H_{0}^{1}(Ω)}: v ∈ S})

Theorem 4.1 says that any subsequence of {u^{k}} has a subsequence which converges strongly in H_{0}^{1}(Ω) to an element of S. In numerical experiments, the entire sequence {u^{k}} converges. Whether or not this is the case is a deep question but we have the following result which indicates that generically there is only one limit point.

Proposition 4.2 Let S(u^{0}) be the set of limit points of Algorithm A_{n} for a given u^{0}. Then either

(i) S(u^{0}) contains exactly one point

or

(ii) S(u^{0}) contains infinitely many points, none of which are isolated.

The key point in the proof of Theorem 4.1 is the following apriori estimate, of independent interest.

Lemma 4.3 Let u ∈ H_{0}^{1}(Ω) be a solution of

\{ \begin{matrix}−Δu = ∑_{i=1}^{n} λ_{i}f′_{i}(u) in Ω \\ F_{i}(u) = ∫_{Ω} f_{i}(u) ≥ γ_{i} > 0, i = 1, …, n \end{matrix} (4.2)

with E(u) ≤ E_{0} < ∞. Then

max_{i} |λ_{i}| ≤ C = C(Ω, E_{0}, γ_{1}^{−1}, …, γ_{n}^{−1}, (Δσ^{−1}))

where Δσ = min_{i} Δσ_{i}.

Assuming the truth of Lemma 4.3 for the moment, let us give the Proof of Theorem 4. As we indicated in the previous section, Lemmas 2.3 and 2.4 permit us to conclude that any subsequence of {u^{k}} has a further subsequence which converges weakly in H_{0}^{1}(Ω) and strongly in L^{2}(Ω) to a solution of (4.2) with

λ_{i} = μ_{i} − τ

Applying Lemma 4.3 we can conclude that λ_{i} > −C, or μ_{i} > τ − C. Hence if we apriori choose a τ > C, we insure that the multipliers μ_{i} > 0, i = 1, … n. This implies, by Lemma 2.4 that F_{i}(u) = γ_{i}, that is, u is a critical point of our original problem P_{n}.

In order to complete the proof of Theorem 4.1, we need to show the strong H_{0}^{1}(Ω) convergence of the subsequence, which we still call u^{k} for convenience. To see this, observe

||u^{k+1} − u||_{H_{0}^{1}(Ω)}^{2} = ∫_{Ω} ∇u^{k+1} · ∇(u^{k+1} − u) − ∫_{Ω} ∇u · ∇(u^{k+1} − u)

= ∑ μ_{i}^{k+1} ∫_{Ω}(u^{k+1} − u)f′_{i}(u^{k}) − τ ∫ u^{k+1}(u^{k+1} − u) − ∫_{Ω} ∇u · (u^{k+1} − u) = ∑ μ_{i}^{k+1} ∫_{Ω}(u^{k+1} − u) f′_{i}(u^{k}) + o(1) as k → ∞ (4.3)

We claim that the μ_{i}^{k+1} are uniformly bounded. For u^{k+1} satisfies

−Δu^{k+1} + τ u^{k+1} = ∑ μ_{i}^{k+1} f′_{i}(u^{k})

and so

\begin{matrix}∑ μ_{i}^{k+1} ∫_{Ω} u^{k}f′_{i}(u^{k}) = ∫ ∇u^{k+1}∇u^{k} + τu^{k+1}u^{k} \\ ≤ E_{τ}(u^{k+1}) + E_{τ}(u^{k} ≤ 2E_{τ}(u^{0}) \end{matrix}

But, ∫_{Ω} u^{k}f′_{i}(u^{k}) ≥ ∫_{Ω} f_{i}(u^{k}) = F_{i}(u^{k}) ≥ γ_{i}, hence ∑_{i} μ_{i}^{k+1}γ_{i} ≤ 2E_{τ}(u^{0}). Since γ_{i} > 0 ∀_{i}, the μ_{i}^{k+1} are uniformly bounded as claimed.

Now observe that for any v ∈ H_{0}^{1}(Ω) ∩ L^{∞}(Ω), v ≥ 0, we have

∑_{i} ∫_{Ω} f′_{i}(v)^{2} dx ≤ λ_{l}^{−1}(Ω) ∑_{i} ∫_{Ω} f′′_{i}(v)|∇v|^{2} = \frac{1}{λ_{1}(Ω)} ∫_{Ω} |∇v|^{2} dx

since f′′_{i}(s) = χ{σ_{i−1} < s < σ_{i}}. Hence,

∑ μ_{i}^{k+1} ∫_{Ω}(u^{k+1} − u)f′_{i}(u^{k}) dx

= ≤ ||u^{k+1} − u||(∑_{i} ∫_{Ω}(μ_{i}^{k+1})^{2}fi^{′2}(u^{k}))^{\frac{1}{2}} dx
≤ c||u^{k+1} − u|| → 0 as k → ∞

(4.4) Therefore, (4.4) (4.5) show that ||u^{k+1} − u||_{H_{0}^{1}(Ω)} → 0 as k → ∞. The proof of Theorem 4.1 is complete. □

We can also give the Proof of Proposition 4.2. Recall from Lemma 2.3 that

E_{τ}(u^{k+1} − u^{k}) ≤ E_{τ}(u^{k}) − E_{τ}(u^{k+1})

and consequently ||u^{k+1} − u^{k}||_{H_{0}^{1}(Ω)} → 0 as k → ∞. Using Theorem 4.1, we can repeat the argument of [], [Proposition 2]. For the readers convenience we will repeat the simple argument. Namely, it suffices to show that if S(u^{0}) contains one isolated point u, then the entire sequence converges to u.

Since u is isolated we can find disjoint neighborhoods N_{1} of u and N_{2} containing all other elements of S(u^{0}). Let 2ϵ = dist(N_{1}, N_{2}). Since ||u^{k+1} − u^{k}||H_{0}^{1}(Ω) ≤ ϵ. We claim that there is an integer k so that for k ≥ k_{1}, u^{k} belongs either to N_{1} or N_{2}. For otherwise, we can find a subsequence u^{k_{i}} lying in the complement of N_{1}∪N_{2}, contradicting Theorem 4.1. Since u ∈ S(u^{0}) we can find a u^{k} ∈ N_{1}, with k ≥ max(k_{0}, k_{1}). The dist(u^{k+1}, N_{2}) ≥ ϵ and so u^{k+1} ∈ N_{1}. By induction u^{j} ∈ N_{1} for j ≥ k and the proposition is proven. □

The remainder of the section will be devoted to the proof of Lemma 4.3. Proof of Lemma 4.3. We first estimate |λ_{1}|, …, |λ_{n−1}|. Set Λ(u) = ∑_{i=1}^{n} λ_{i}f′_{i}(u) and let Λ_{i} = Λ(σ_{i}), i = 0, …, n − 1. Then

λ_{0} = 0 and λ_{i} = Λ_{i−1} + λ_{i}Δσ_{i}Λ(s) = Λ_{i−1} + λ_{i}(s − σ_{i−1}) on [σ_{i−1}, σ_{i}]

= Λ_{i} + λ_{i}(s − σ_{i}) (4.5)

Since λ_{i} = (Λ_{i} − Λ_{i−1})/Δσ_{i}, it suffices to bound |Λ_{1}|, …, |Λ_{n−1}|. We fix the index i so that |Λ_{i}| = max |Λ_{j}|_{1≤j≤n−1} (4.6)

Then for s ∈ I =: [σ_{i} − \frac{1}{3}Δσ_{i}, σ_{i} − \frac{1}{6}Δσ_{i}],

|Λ(s)| ≥ \frac{1}{3}|Λ_{i}| (4.7)

since |Λ(s)| ≥ |λ_{i}| · \frac{1}{3}Δσ_{i} and |λ_{i}| ≤ 2\frac{|Λ_{i}|}{Δσ_{i}}.

Let ϕ(s) be a Lipschitz cut off function with ϕ(s) = 1 for s ∈ I, ϕ(s) = 0 for s ≤ σ̄_{i} = \frac{1}{2}(σ_{i−1} + σ_{i}) and s ≥ σ_{i} and 0 ≤ ϕ(s) ≤ 1. Using ϕ(u) as a test function in (4.21) we have the identity

∫_{Ω} ϕ′(u)|∇u|^{2} dx = ∫_{Ω} ϕ(u)Λ(u) dx (4.8)

Since Λ(s) has one sign on the support of ϕ(s), it follows from (4.9) that

\begin{matrix}∫_{\{u∈I\}} |Λ(u)| dx ≤ ∫_{Ω} ϕ(u)|Λ(u)| dx ≤ ∫_{Ω} |ϕ′(u)| |∇u|^{2} dx \\ ≤ \frac{6}{Δσ_{i}} ∫_{Ω} |∇u|^{2} dx = 12(Δσ_{i})^{−1}E(u) \end{matrix} (4.9)

Combining (4.10) and (4.8), we have

|Λ_{i}| ≤ \frac{36E(u)}{(Δσ_{i})|\{u ∈ I\}|} (4.10)

Thus to estimate |Λ_{i}| we need a positive lower bound for |{u ∈ I}|.

Let σ = σ_{i} − \frac{1}{3}Δσ_{i}, σ′ = σ_{i} − \frac{1}{6}(Δσ)_{i} so that I = [σ, σ′] ⊂ (σ_{i−1}, σ_{i}), |I| = \frac{1}{6}Δσ_{i} for some 1 ≤ i ⊆ n − 1. Define

f_{I}(s) = \frac{1}{2}(s − σ)_{+}^{2} − (s − σ′)_{+}^{2})

Then

||f′_{I}(u)||L^{\frac{n}{n−1}_{Ω}} ≥ |I| |\{u > Γ_{n−1}\}|^{\frac{n−1}{n}} (4.11)

On the other hand, by the Sobolev inequality

||f′_{I}(u)||L^{\frac{n}{n−1}_{(Ω)}} ≤ c ∫_{Ω} f′′_{I}(u)|∇u| dx ≤ c|\{u ∈ I\}|E(u)^{\frac{1}{2}} (4.12)

and by the Poincaré inequality. 2F_{n}(u) = ||(u − σ_{n−1})_{+}||_{L^{2}(Ω)}^{2} ≤ C_{n}|\{u > σ_{n−1}\}|^{\frac{2}{N}} E(u). (4.13) Combining (4.12) - (4.14) gives the required lower bound |{u ∈ I}| ≥ c_{N}|I|^{2} F_{n}(u)^{N−1} E(u)^{−N}. (4.14) Together with (4.11) this gives the estimate |Λ_{i}| ≤ C_{N}E(u)^{N+1}F_{n}(u)^{1−N}(Δσ_{i})^{−3} ≤ C_{N}E_{0}^{N+1}γ_{n}^{1−N}(Δσ)^{−3}. (4.15)

To complete the proof we need only bound |λ_{n}|. Consider −Δu = Λ_{n−1} + λ_{n}(u − σ_{n−1}) in {u > σ_{n−1}} (4.16) Using (u − σ_{n−1})_{+} as a test function in (4.17) yields the identity ∫_{Ω} |∇(u − σ_{n−1})_{+} + 2λ_{n}F_{n}(u) (4.17) Therefore,

γ_{n}|λ_{n}| ≤ E_{0} + |Λ_{n−1}| · (\frac{|Ω|}{λ_{1}(Ω)}E_{0})^{\frac{l}{2}} (4.18)

Recalling (4.16) this gives |λ_{n}| ≤ C(Ω, N_{1}E_{0})γ_{n}^{−2}(Δσ)^{−3} (4.19) This completes the proof of Lemma 4.3. □

5 Generalizations.

In this section we give three separate extensions of the prototype problems examined in the preceding sections. Since each extension is fairly straightforward our discussion here will be brief.

First, we introduce a free-boundary into problems (P_{∞}) and (P_{n}). This is accomplished simply by informing the constraints (2.4c) only on the interval σ_{0} ≤ σ < + ∞ for some given σ_{0} > 0. Thus, the problem (P_{∞}) becomes

(5.1) E(u) → min subject to ∫_{Ω}(u − σ)_{+} dx = β(σ), σ_{0} ≤ σ < + ∞.

Similarly, the problem (P_{n}) is defined with respect to a partition 0 < σ_{0} < σ_{1} < …, < σ_{n−1} < σ_{n} = + ∞, and the constraints for (P_{n}) are defined by (2.6), (2.7) as before. Now, however, f_{i}(s) = 0 for all s ≤ σ_{0} and i = 1, …, n. the validity of Theorem 2.2 remains unaltered in this extension, and a solution u of (P_{n}) satisfies −Δu = Λ(u) in Ω, u = 0 on ∂Ω, where

Λ(s) = \{ \begin{matrix}0 & if s < σ_{0} \\ ∑_{i=1}^{n} λ_{i} f′_{i}(s) & if s ≥ σ_{0} \end{matrix}

This may be interpreted as the equivalent free-boundary problem

−Δu = \{ \begin{matrix}0 & in \{0 < u < σ_{0}\} \\ Λ(u) & in \{u > σ_{0}\} \end{matrix} (5.2)

|∇u| is continuous across {u = σ_{0}},

the latter being the free-boundary condition. As a simple illustration, we note that the case of one constraint (n = 1) yields

\{ \begin{matrix}−Δu = λ_{1}(u − σ_{0})_{+} in Ω, u = 0 on ∂Ω \\ ∫_{Ω} (u − σ_{0})_{+}^{2} = 2γ_{1}, \end{matrix}

which is a familiar free-boundary problem.

The algorithm (A_{n}) constructed in §3 and its convergence properties established in §4 admit immediate generalizations in the extension just described. Indeed, no essential changes in the above development are needed to include the free-boundary into (P_{n}).

Next, we consider an extension of (P_{∞}) and (P_{n}) (without free-boundaries) which allows u and ū to change sign in Ω. The requirement that ū ≥ 0 in Ω is now relaxed, and the constraints in (P_{∞}) are replaced by

\{ \begin{matrix}∫_{Ω}(u − σ)_{+} dx = β^{+}(σ) for 0 ≤ σ < + ∞ \\ ∫_{Ω}(u − σ)_{−} dx = β^{−}(σ) for − ∞ < σ ≤ 0, \end{matrix}

where β^{+}, β^{−} are defined by ū (so that the above holds with ū substituted for u). The extension of (P_{n}) is then defined with respect to two partitions 0 = σ_{0}^{+} < σ_{1}^{+} < … < σ_{n−1}^{+} < σ_{n}^{+} = + ∞ and 0 = σ_{0}^{−} > σ_{1}^{−} > … > σ_{m−1}^{−} > σ_{m}^{−} = − ∞ for some n and m. The functionals F_{i}^{+}(i = 1, …, n) and F_{j}^{−}(j = 1, …, m) corresponding to these partitions are defined as before by the convex functions

f_{i}^{+}(s) = \frac{1}{2}(s − σ_{i=1}^{+})_{+}^{2} − \frac{1}{2}(s − σ_{i}^{+})_{+}^{2}

f_{j}^{−}(s) = \frac{1}{2}(s − σ_{j−1}^{−})_{+}^{2} − \frac{1}{2}(s − σ_{j}^{−})_{+}^{2};

also, the constraint values are taken to be

γ_{i}^{+} = ∫_{σ_{i−1}^{+}}^{σ_{i}^{−}} β^{+}(σ) dσ, γ_{j}^{−} = ∫_{σ_{j}^{−}}^{σ_{j−1}^{−}} β^{−}(σ) dσ.

The constraints for (P_{n}) are therefore replaced by the n + m constraints F_{i}^{+}(u) = γ_{i}^{+}, F_{j}^{−}(u) = γ_{j}^{−}, and so an analogous variational problem is defined, which we shall call (P_{n,m}). If each γ_{i}^{+} > 0 and each γ_{j}^{−} > 0, then the analogue of Theorem 2.2 holds, and a solution of (P_{n,m}) satisfies −Δu = Λ(u) in Ω, u = 0 on ∂Ω, where now

Λ(s) = ∑_{i=1}^{n} λ_{i}^{+} f_{i}^{+′}(s) + ∑_{j=1}^{m} λ_{j}^{−} f^{−′}(s),

a piecewise linear function on −∞ < s < + ∞ which is increasing whenever λ_{1}^{+} ≥ 0 and λ_{j}^{−} ≥ 0 for every i and j.

Again the algorithm and its convergence properties require no essential changes to handle the problem (P_{n,m}).

Finally, we remark that a variable coefficient version of (P_{n}) can be treated as a straightforward extension of the prototype problem. In this version the objective and constraint functionals are taken to be

E(u) = \frac{1}{2} ∫_{Ω}\{ ∑_{p, q=1}^{N} a_{pq}(x)u_{x_{p}} u_{x_{q}} + a(x)u^{2}\} dx

F_{i}(u) = ∫_{Ω} b(x)f_{i}(u) dx,

where the functions a_{pq}(x), a(x), b(x) belong to C^{α}(Ω̄)(0 < α < 1) and satisfy the conditions

∑_{p,q=1}^{N} a_{pq}(x)ξ_{p}ξ_{q} ≥ θ|ξ|^{2} for all ξ ∈ ^{N}

a(x) ≥ 0, b(x) ≥ θ′

uniformly for x ∈ Ω for some positive constraints θ and θ′. Then a solution of (P_{n}) solves the nonlinear elliptic eigenvalue problem.

\{ \begin{matrix}− ∑_{p,q} \frac{∂}{∂x_{p}} (a_{pq}(x)\frac{∂u}{∂x_{q}}) + a(x)u = b(x)λ(u) in Ω \\ u = 0 on ∂Ω, \end{matrix}

where Λ(s) is as before.

It is straightforward to check that our development in §2-4 as well as the two preceded extensions can be carried out in this more general context.

6 Appendix 1. Variational Problems in Magneto Hydrodynamics.

Equilibrium problems in magneto hydrodynamics (MHD) supply some of the main examples of multiconstrained variational problems of the kind that we study in this paper. These model problems in plasma physics arise both in controlled thermonuclear fusion research and in astrophysics. The classical approach to solving these problems makes use of a general variational principle which characterizes equilibrium configurations of a plasma assumed to be governed by ideal MHD. This principle depends upon determining the complete family of quantities (expressed as volume integrals) which are conserved under the governing evolution equations. Then it identifies an equilibrium as a minimizer of the total energy over the class of configurations which maintain given values of all of the other conserved quantities. Although this general characterization can be stated formally for fully three-dimensional configurations, we shall consider it only under an assumption of spatial symmetry - either two-dimensionality, axial symmetry or helical symmetry - since then it takes a simpler and more tractable form. This particular form of the variational principle for symmetric equilibrium configurations has been given by Woltjer [?, ?], who has also determined the complete family of conserved quantities. (He emphasizes the axially symmetric case, but his analysis can be modified to apply as well to the other symmetric cases.) We shall briefly describe the case of magneto static equilibrium (for which there is no mass flow) in two-dimensions; we shall leave aside the minor technical modifications needed to treat axial or helical symmetry since, even though these cases are important in real applications, they are the same conceptually.

We begin by posing an abstract variational problem similar to (P_{∞}) which encompasses the physical problems of interest as special cases. Let Ω ⊆ ^{2} be the cross-section of a cylindrical domain Ω × , and let (x_{1}, x_{2}) denote the variable point in Ω and x_{3} be the ignorable coordinate. Consider the minimization problem

∫_{Ω} [\frac{1}{2}|∇u|^{2} + h_{1}(v_{1}), +h_{2}(v_{2})] dx → min over
∫_{Ω} v_{1}(u − σ)_{+} dx = β_{1}(σ)
∫_{Ω} v_{2}(u − σ)_{+} dx = β_{2}(σ), (σ_{0} ≤ σ < + ∞)

(6.1)

where the admissible triple (u, v_{1}, v_{2}) belongs to H_{0}^{1}(Ω) × L^{r_{1}}(Ω) × L^{r_{2}}(Ω) for some 1 < r_{1}, r_{2} < + ∞. The given functions h_{1} and h_{2} are assumed to be smooth and strictly convex with h_{ℓ}(0) = h′_{ℓ}(v) = 0 and h_{ℓ}(z) = O(|z|^{r_{ℓ}}) as |z| → ∞ (ℓ = 1, 2). As in (P_{∞}), the two infinite family of constraints are parametrized by σ ∈ [σ_{0}, + ∞), and β_{1}(σ) and β_{2}(σ) are given data.

The physical interpretation of (6.1) is as follows. The magnetic field \underline{B} = (B′, B^{2}, B^{3}), which is independent of x_{3}, satisfies ∇ · B = 0 in Ω × and hence admits a representation \underline{B} = (u_{x_{2}}, −u_{x_{1}}, v_{1}), where u is the flux function (or stream function) for its poloidal part and v_{1} is its toroidal part. The magnetic energy density (per unit volume) is then \frac{1}{2}|\underline{B}|^{2} = \frac{1}{2}|∇u|^{2} + \frac{1}{2}v_{1}^{2}. The mass density ρ of the plasma is represented by v_{2}. The internal energy density (per unit volume) is given by ρ^{γ}/(γ − 1), in accordance with the polytropic law p = ρ^{γ} with p denoting pressure. Therefore, the objective functional represents total (potential) energy when we put

h_{1}(v_{1}) = \frac{1}{2}v_{1}^{2}, h_{2}(v_{2}) = \frac{v_{2}^{γ}}{γ − 1} (6.2)

The interpretation of the constraints relies on differentiating them with respect to the parameter σ, then there results

∫_{{u>σ}} v_{1} dx = −β′_{1}(σ), ∫_{{u>σ}} v_{2} dx = −β′_{2}(σ)

All of these integrals are extended over the interior of a (cylindrical) flux surface {u = σ} - that is, a flux tube {u > σ}. But ideal MHD requires that, in evolution, each flux tube must move with the flow preserving its strength and mass. It is readily verified that the above integrals are, respectively, the toroidal flux and mass of the flux tube {u > σ}, and hence they are conserved quantities. (The conservation of poloidal fluxes is implicit in the parametrization which uses the values of the flux function u.)

As is shown in Lemma 2.1, the constraints in (6.1) imply that all constraints involving integrals of the form

∫_{Ω} v_{ℓ}ϕ(u) dx where ϕ is an arbitrary function.

Indeed, the latter integrals are used to express the family of conserved quantities in [?, ?]. Roughly speaking, the functions (s − σ)_{+} (σ_{0} ≤ σ < + ∞) are chosen here as a particularly useful basis for the space of all functions ϕ(s).

This allows us to give a more direct physical interpretation to these integrals. Moreover, this also permits us to discretize the constraints in precisely the same way as the prototype problem (P_{n}) is formed from (P_{∞}). Recalling the definition (2.7), we replace the infinite family of constraints in (6.1) by

\{ \begin{matrix}∫_{Ω} v_{1} f_{i}(u) dx = γ_{1i} \\ ∫_{Ω} v_{2} f_{i}(u) dx = γ_{2i}, (i = 1, …, n) \end{matrix}

relative to a (fixed) partition 0 ≤ σ_{0} < σ_{1} < … < σ_{n−1} < σ_{n} = + ∞. The resulting discretized version of the minimization problem (6.1) is now accessible to analysis and computation. Its variational equations can be calculated formally to be

\{ \begin{matrix}−Δu & = & v_{1}Φ′1(u) + v_{2}Φ′_{2}(u) \\ h′_{1}(v_{1}) & = & Φ_{1}(u) \\ h′_{2}(v_{2}) & = & Φ_{2}(u), \end{matrix} (6.3)

where, by virtue of the Lagrange multiplier rule,

Φ_{ℓ}(s): = ∑_{i=1}^{n} λ_{ℓi} f_{i}(s) (ℓ = 1, 2) (6.4)

for some multiplier λ_{ℓi}. The observation that v_{1} and v_{2} enter into these equilibrium equations algebraically suggests immediately that they be eliminated. This yields an equation for u alone:

−Δu = (h′_{1}(Φ_{1}(u))Φ′_{1}(u) + (h′_{2})^{−1}(Φ_{2}(u))Φ′_{2}(u), (6.5)

which is a generalized form of the familiar Grad-Shafranov equation. In order to express this equation more plainly, it is useful to introduce the convex conjugate functions to h_{ℓ}(ℓ = 1, 2) defined by h_{ℓ}(z^{∗}) = sup_{z} zz^{∗} − h_{ℓ}(z). Then (6.5) can be written as simply

−Δu = P′_{1}(u) + P′_{2}(u) with P_{ℓ}(s): = h_{ℓ}^{∗} (Φ_{ℓ}(s)). (6.6)

Returning to the physical case, we find that

h_{1}^{∗}(z^{∗}) = \frac{1}{2}(z^{∗})^{2}, h_{2}^{∗}(z^{∗}) = (\frac{γ − 1}{γ}z^{∗})^{\frac{γ}{γ−1}} ,

and consequently we obtain the relations

P_{1} = \frac{1}{2}v_{1}^{2}, P_{2} = v_{2}^{γ}.

It follows from these relations that the profile functions P′_{1} and P′_{2} in (6.6) have the usual interpretations or current and pressure profiles, respectively, in the Grad-Shafranov equation.

Of course, the principal novelty inherent in (6.6) lies in the fact that the functions P′_{ℓ} are not specified, but are instead determined along with the solution u through (6.4). We see therefore that (6.6) is a “generalized (or queer) differential equation” (GDE) in the sense of Grad et.al.[?, ?]. However, our viewpoint differs from that proposed by Grad, who viewed a GDE as a highly implicit equation for u which combined the elliptic operation Δ and differentiation with respect to the volume variable α = α(σ) = |{u > σ}|. Rather, we recognize (6.6) as merely the v_{1}, v_{2} - eliminated form of the standard variational equations for a solution triple (u, v_{1}, v_{2}) of the classical minimization problem (6.1). Thus, by treating that variational problem directly we obviate the need to introduce the notion of a GDE. On the other hand, the compelling reasons put forward by Grad and others for prescribing conserved quantities associated with the evolution equations governing ideal MHD in place of the current and pressure profiles in the Grad-Shafranov equation retain their validity; indeed, they provide the main justification for our multi-constrained variational approach.

The relevant physical problem also involves a free-boundary, the interface between the (confined) plasma and the (surrounding) vacuum. This feature can be introduced (as in §5) by fixing σ_{0} > 0 in (6.1). Then the free-boundary is {u = σ_{0}}, and {u > σ_{0}} and {0 < u < σ_{0}} are the plasma and vacuum regions, respectively. The flux constant σ_{0} effectively specifies the (total) poloidal flux in the vacuum region. an additional constraint ∫_{Ω} v_{1} dx = γ_{10} can be imposed to specify the (total) toroidal flux in the vacuum region. (This is necessary since the basis functions f_{i}(s) vanish for s ≤ σ_{0}). The mass density v_{2} will be positive in {u > σ_{0}} and vanish in {0 < u < σ_{0}} whenever the solution satisfies λ_{2_{i}} ≥ 0 (i = 1, …, n), and this condition may be expected to hold in general. The free-boundary, ∂{v_{2} > 0}, is therefore a flux surface on which pressure vanishes, as required by the physical interface conditons.

We conclude this discussion by noting some special cases of the above development. First, we consider the case of an incompressible fluid (or plasma) which we obtain by setting v_{2} = 1 (uniform density) in (6.1). Equivalently, we may achieve this limit case by taking the exponent γ to infinity; then we can verify that the compressible solutions tend to an incompressible solution in the limit. Second, we have the case of a purely poloidal magnetic field which we get by putting v_{1} = 0. This specialization of (6.1) results in a substantial simplification since one family of constraints is dropped. Third, we combine the two preceding cases and we arrive at the prototype problem (P_{∞}) or its discretized form (P_{n}), both posed in terms of u alone.

With this result in mind we remark that it is possible to construct algorithms analogous to (A_{n}) which iteratively solve the more general problem (6.1), at last under favorable circumstances. However, a convergence theory similar to that given in §3-5 for the prototype problem is not yet available for this general problem, owing mainly to the lack of joint convexity of the constraint functionals in (6.1).

7 Appendix 2. Some Comments On Algorithm (A_{n}).

This appendix is included in order to indicate the main ideas underlying the algorithm (A_{n}) presented formally in the preceding section. As such the comments given here tend to emphasize the conceptual rather than logical development of general iterative procedure for solving the variational problem (P_{n}).

The globally convergent iterative procedure (A_{n}) can be considered as the product of a development consisting of three stages. At each of these stages a key ingredient in the algorithm is introduced. Therefore, it is useful when explaining the algorithm (A_{n}) to describe the form that it takes at each stage of its development. We now proceed to give such a description, stressing the ideas involved rather than the details of proof.

The first stage of (A_{n}) consists, in essence, of an iterative procedure which is based upon linearizing the nonlinear equality constraints occurring in the minimization problem (P_{n}). This first form of (A_{n}) is defined iteratively by the equations

−Δu^{k+1} = ∑_{j=1}^{n} λ_{j}^{k+1}f′_{j}(u^{k}) in Ω, u^{k+1} = 0 on ∂Ω (7.1)

∫_{Ω} [f_{i}(u^{k}) + f′_{i}(u^{k})(u^{k+1} − u^{k})] dx = γ_{i} (i = 1, …, n). (7.2)

For a fixed initialization u^{0} ∈ M_{n}(ū), this iterative procedure defines a sequence (u^{k+1}, λ^{k+1}) ∈ H^{1}0(Ω) × ^{n} which is intended to converge (in some suitable sense) to a solution of (P_{n}). Moreover, since the defining equations (7.1) and (7.2) are linear in (u^{k+1}, λ^{k+1}), the explicit construction of that iterative pair from u^{k} is a routine matter. First, we solve the linear elliptic problems

−Δw_{j} = f′_{j}(u^{k}) in Ω, w_{j} = 0 on ∂Ω (j = 1, …, n); (7.3)

second, we find λ^{k+1} ∈ ^{n} by solving the n × n system of linear algebraic equations ∑ a_{ij}(u^{k})λ_{j}^{k+1} = c_{i}(u^{k}) (i = 1, …, n), (7.4) where, as in (2.29) and (2.30) with τ = 0,

a_{ij}(u^{k}) = ∫_{Ω} ∇w_{i} · ∇w_{j} dx c_{i}(u^{k}) = γ_{i} − ∫_{Ω} [f_{i}(u^{k}) − f′_{i}(u^{k})u^{k}] dx;

third, we set

u^{k+1} = ∑_{j=1}^{n} λ_{j}^{k+1} w_{j}. (7.5)

This algorithm treats the constraints in (P_{n}) which the main source of difficulty in the problem, by replacing them at each iterative step with their linearization about the previous iterate. If the iterative sequence (u^{k+1}, λ^{k+1}) converges to (u, λ) ∈ H_{0}^{1}(Ω) ×^{n} (in norm, say ), then u is a critical point of (P_{n}), or equivalently, (u, λ) solves (S), meaning that

∫_{Ω} ∇u · ∇v dx = ∑_{j=1}^{n} λ_{j} ∫_{Ω} f′_{j}(u)v dx (7.6)

∫_{Ω} f_{i}(u) dx = γ_{i} (i = 1, …, n) (7.7)

for all v ∈ H_{0}^{1}(Ω). Thus, the constraints (7.7) on solutions are retrieved from the linearized constraints (7.2) on iterates. Furthermore, it is noteworthy that the defining equations (7.1), (7.2) are precisely the variational equations for the convex minimization problem

E(u) → min subject to

F_{i}(u^{k}) + ⟨F′_{i}(u^{k}), u − u^{k}⟩ = γ_{i} (i = 1, …, n); (7.8)

and it is obvious that u^{k+1} is the unique minimizer for (7.8) and λ^{k+1} is the uniquely determined multiplier vector for u^{k+1}.

Although the problem (P_{n}) is itself nonconvex having constraints that are nonlinear equalities, the convexity of its objective and constraint functionals implies that under certain conditions the sequence of iterates enjoys special monotonicity and convergence properties. One consequence of these convexity attributes is the inequality.

∫_{Ω} f_{i}(u^{k+1}) dx ≥ ∫_{Ω} [f_{i}(u^{k}) + f′_{i}(u^{k} − u^{k})] dx = γ_{i} (7.9)

which holds unconditionally for every k and each i = 1, …, n. Hence, the iterate u^{k+1} may be said to satisfy the constraints (7.7) in a relaxed sense with inequalities replacing the equalities. From (7.9) combined with Lemma 4.2 it follows that the sets {x ∈ Ω: σ_{i−1} < u(x) < σ_{i}} have positive measure. Using the argument given in the proof of Theorem 2.2 the positive-definiteness of the n × n matrix A = (a_{ij}(u^{k})) is therefore ensured, since for any ξ ∈ ^{n}

∑_{i, j=1}^{n} a_{ij}(u^{k})ξ_{i}ξ_{j} = ∫_{Ω} |∇(ξ · w|^{2} dx > 0

unless ξ · w = 0 in Ω, which is equivalent to ξ · f′(u^{k}) = 0 in Ω. Thus, the invertibility of A is derived from (7.9), and in turn the explicit construction of the iterate (u^{k+1}, λ^{k+1}) is seen to be well-defined. Another consequence of convexity is the monotonicity of the values of the objective functional along the iterative sequence, but, unlike (7.9), this property requires that a further condition be satisfied–namely, λ_{i}^{k+1} ≥ 0 for each i = 1, …, n. A simple calculation making use of (7.1), (7.2) and (7.9) yields the result: E(u^{k}) − E(u^{k+1}) − E(u^{k+1} − u^{k}) = ∫_{Ω} ∇u^{k+1} · ∇(u^{k} − u^{k+1}) dx

= ∑_{j=1}^{n} λ_{j}^{k+1} ∫_{Ω} f′_{j}(u^{k})(u^{k} − u^{k+1}) dx

= ∑_{j=1}^{n} λ_{j}^{k+1}{∫_{Ω} f_{j}(u^{k}) dx − γ_{j}}
≥ 0.

Therefore,

E(u^{k+1}) ≤ E(u^{k}) and E(u^{k+1} − u^{k}) → 0 as k → + ∞, (7.10)

provided that every λ_{i}^{k+1} is nonnegative. Moreover, this condition also furnishes an upper bound on |λ^{k+1}|, as is evident from the estimate

E(u^{k+1}) + E(u^{k}) ≥ ∫_{Ω} ∇u^{k+1} · ∇u^{k} dx = ∑_{j=1}^{n} λ_{j}^{k+1} ∫_{Ω} f′_{j}(u^{k})u^{k} dx

≥ ∑_{j=1}^{n} λ_{j}^{k+1} γ_{j}.

These properties lead directly to a (generalized) convergence result for the iterative sequence. Specifically, it follows that every subsequence of (u^{k+1}, λ^{k+1}) has a further subsequence which converges to a solution u, λ) of (7.6), (7.7), with the subsequence of u^{k+1} converging to u weakly in H_{0}^{1} and strongly in L^{2}. In other words, the iterative procedure produces a sequence which converges (in the sense of minimum distance) to the set S of solutions whenever the condition λ_{j}^{k+1} ≥ 0 can be enforced.

One would expect that the latter condition could be satisfied if the initialization u^{0} of the above iterative procedure is taken close enough to a solution of (P_{n}) which itself satisfies λ_{i} ≥ δ > 0(i = 1, …, n) for some δ. In addition, the monotonicity property (7.10) suggests that, at least typically, the iterates should converge to a local minimum of E rather then merely a solution of S. Thus, the utility of this provisional form of the algorithm for computing solutions of (P_{n}) is demonstrated within certain limits. On the other hand, the condition that λ_{j}^{k+1} ≥ 0 is quite stringent and cannot be expected necessarily to hold when either (i) u^{0} is not close enough to a local minimizer or (ii) a solution u for which some λ_{i} are negative is sought.

The second stage of (A_{n}) is a modification of the above iterative procedure which weakens the above condition. Now the parameter τ > 0 is introduced and (u^{k+1}, λ^{k+1}) is defined by the analogues of equations (7.1), (7.2) where −Δ is replaced by −Δ+τ and λ_{j}^{k+1} is replaced by λ_{j}^{k+1}+τ. With these substitutions made throughout the above algorithm, it can be verified that the analogues of the above mentioned monotonicity and convergence properties hold provided that λ_{i}^{k+1} ≥ −τ for each i = 1, …, n; with these modification it can be seen that E is replaced by E_{τ} = E + τ(F_{1} + ⋯ + F_{n}), which is exactly the functionals defined in (2.20). This device for shifting the iterative multipliers by τ rests crucially on the identity (2.14), and hence on the structure of the family of constraints imposed in (P_{n}). Presumably, with this second form of (A_{n}) in hand one would expect to be able to treat the general case of (P_{n}) provided that the initialization of the iterative procedure is close enough to a solution. The choice of a suitable τ would then be dictated by the magnitude of the negative λ_{i} values corresponding to the solution.

The third stage for (A_{n}) is precisely the content of §3. The final form of the algorithm differs from the preceding forms in that the (linearized) constraints imposed at each iteration are now relaxed to inequalities. As is shown in Lemma 2.3 and 2.4, the monotonicity and (partial) convergence properties of this iterative procedure hold unconditionally, since μ_{i}^{k+1} ≥ 0 by necessity. Therefore, the algorithm (A_{n}) defines a sequence (u^{k+1}, λ^{k+1}) which converges in the sense explained above given as arbitrary initialization u^{0}. This means that (A_{n}) is a globally convergent iterative procedure, in contrast to its two provisional forms with equality constraints. In essence, the reason for this property of (A_{n}) lies in the fact that it takes full advantage of the convexity attributes of the objective and constraint functionals.

It is necessary however to show that the resulting solution u satisfies the constraints (7.7) as equalities. This is accomplished by choosing τ large enough depending on the solutions of S (not on the iterates) so that λ_{i} + τ > 0 for each i, and hence that each constraint must be active.