Numerical Maximization of the p-Laplacian Energy of a Two-Phase Material

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. SIAM J. NUMER. ANAL.©2021 Society for Industrial and Applied Mathematics Vol. 59, No. 6, pp. 3077--3097 NUMERICAL MAXIMIZATION OF THE p-LAPLACIAN ENERGY OF A TWO-PHASE MATERIAL\ast JUAN CASADO-D\' IAZ\dagger , CARLOS CONCA\ddagger ,AND DONATO V\' ASQUEZ-VARAS\S Abstract. For a diffusion problem modeled by the p-Laplacian operator, we are interested in obtaining numerically the two-phase material which maximizes the internal energy. We assume that the amount of the best material is limited. In the framework of a relaxed formulation, we present two algorithms, a feasible directions method and an alternating minimization method. We show the convergence for both of them, and we provide an estimate for the error. Since for p > 2 both methods are only well-defined for a finite-dimensional approximation, we also study the difference between solving the finite-dimensional and the infinite-dimensional problems. Although the error bounds for both methods are similar, numerical experiments show that the alternating minimization method works better than the feasible directions one. Key words. optimal design, two-phase material, p-Laplacian operator, feasible directions method, alternating minimization method AMS subject classifications. 49M05, 49J20 DOI. 10.1137/20M1353563 1. Introduction. The aim of the present work is the numerical resolution of an optimal design problem. It corresponds to the maximization of the energy for a nonlinear diffusion process in a two-phase material modeled by the p-Laplacian operator. Namely, we are interested in the control problem (1.1) \left\{ max \omega 1 p\int \Omega (\alpha \scrX \omega +\beta (1 - \scrX \omega )) | \nabla u| pdx - div \bigl( (\alpha \scrX \omega +\beta (1 - \scrX \omega )) | \nabla u| p - 2\nabla u\bigr) =fin \Omega u\in W1,p 0(\Omega ), \omega \subset \Omega measurable, | \omega | \leqslant \kappa , with \Omega a bounded open set in \BbbR N,N\geqslant 2, p\in (1,\infty ), \alpha , \beta , \kappa > 0, \alpha < \beta , and f\in W - 1,p\prime (\Omega ). Here \alpha and \beta are the diffusion constants corresponding to the two materials that we want to mix in order to maximize the corresponding functional. If we do not impose any restrictions on the amount of material \alpha (i.e., \kappa \geqslant | \Omega | ), then the solution is the trivial one given by \omega = \Omega . Thus, the interesting case corresponds to \kappa < | \Omega | . This problem has been extensively studied for p= 2 ([1], [6], [8], [12], \ast Received by the editors July 17, 2020; accepted for publication (in revised form) August 23, 2021; published electronically December 21, 2021. https://doi.org/10.1137/20M1353563 Funding: The first author has been partially supported by Project MTM2017-83583 of the Ministerio de Ciencia, Innovaci\'on y Universidades of Spain. The second author is partially supported by PFBasal-001 (CeBiB), FB210005 (CMM), and ACE210010, BASAL funds for centers of excellence from ANID-Chile and by Regional Program NEMBICA STIC190013 STIC-AmSud. The third author has been partially supported by the CONICYT PFCHA/DOCTORADO BECAS CHILE/2018-21182101. \dagger Departamento de Ecuaciones Diferenciales y An\'alisis Num\'erico, Universidad de Sevilla, Sevilla 41012, Spain (jcasado[email protected]). \ddagger Department of Engineering Mathematics, Center for Mathematical Modelling (CMM), UMI 2807 CNRS-Chile, and Center for Biotechnology and Bioengineering (CeBiB), University of Chile, Santiago, RM, 837 0459, Chile ([email protected]). \S Center for Mathematical Modelling (CMM), UMI 2807 CNRS-Chile, University of Chile, Santiago, RM, 837 0459, Chile (dv[email protected]). 3077 Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 3078 J. CASADO-D\' IAZ, C. CONCA, AND D. V\' ASQUEZ-VARAS [14], [15], [20]). In this case, it models, for example, the optimal rearrangement of two materials in the cross section of a beam in order to minimize its torsion (in this application f= 1). Analogously, for p\in (1,2) \cup (2,\infty ), the p-Laplacian operator models the torsional creep in the cross section of a beam [9]. Therefore, problem (1.1) corresponds to find the two-phase material which minimizes the torsion in nonlinear elasticity, assuming that the amount of the best material is limited. As usual for this type of problem ([17], [18]), it has no solution in general ([6], [7], [20]). Thus, it is necessary to work with a relaxed formulation which can be obtained from the homogenization theory ([1], [19], [21]). In the present case, it has been proved in [7] ([20] for p= 2) that such relaxation is given by (1.2) \left\{ max u,\theta 1 p\int \Omega | \nabla u| p (1 + c\theta )p - 1dx - div\Bigl( | \nabla u| p - 2 (1 + c\theta )p - 1\nabla u\Bigr) =1 \beta fin \Omega u\in W1,p 0(\Omega ), \theta \in L\infty (\Omega ; [0,1]),\int \Omega \theta dx \leqslant \kappa , with c= (\beta \alpha )1 p - 1 - 1. In this formulation, the materials \alpha and \beta have been replaced by mixtures of them obtained by laminations. The new control variable \theta represents the proportion of the best material \alpha used in the mixture. The problem can also be formulated in a simple way as the following calculus of variation problem: (1.3) \left\{ min u,\theta \biggl\{ 1 p\int \Omega | \nabla u| p (1 + c\theta )p - 1dx - 1 \beta \langle f, u\rangle \biggr\} \theta \in L\infty (\Omega ; [0,1]), u \in W1,p 0(\Omega ),\int \Omega \theta dx \leqslant \kappa . The numerical resolution of (1.3) for p= 2 has been the subject of several works. In this way, some numerical simulations have been carried out in [12] and [14] using a multigrid method. In [1] and [22], the convergence of the alternating minimization algorithm has been shown, using the optimality conditions. In [3], the convergence of a projected gradient method has been studied. For p\not = 2, the use of the optimality conditions implies the resolution of the p- Laplacian equation in each iteration. This is a problem which has been considered, for example, in [11] and [13] using a steepest descent method. We also refer the reader to [10], where a reformulation of the p-Laplacian is given in order to use an augmented Lagrangian method. In these works, the order of convergence is linear in the best case. In the present paper, we introduce two algorithms to solve (2.1). The first one is based on the Frank--Wolfe algorithm, also known as the feasible direction method. The second one is an alternating minimization method. In both of them, we choose a descent direction in H1 0(\Omega ) instead of W1,p 0(\Omega ) and we solve a linear problem instead of a p-Laplacian, which, as we said above, is very expensive from a computational point of view. For p > 2, this forces us to work with a discretized version of the problem because H1 0(\Omega ) is not contained in W1,p 0(\Omega ). We prove the convergence of both methods obtaining estimates for the rate of convergence. In the best of the cases (p\geqslant 2), we only have a convergence of order 1/i, with ithe number of iterations. This is due to the nonstrict convexity of the problem. In this sense, we can observe that solving the minimum in \theta in problem Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. NUMERICAL MAXIMIZATION OF THE p-LAPLACIAN ENERGY 3079 (1.3) and using the Kuhn--Tucher theorem, we can rewrite (1.3) as (see [6], [7], [12], [14]) (1.4) max \mu \geqslant 0min u\in W1,p 0(\Omega ) \biggl\{ \int \Omega F(\mu , \nabla u)dx - \mu p(p - 1)c p\kappa - 1 \beta \langle f, u\rangle \biggr\} with F(\mu , \xi ) = \left\{ 1 p | \xi | p (1 + c)p - 1+\mu p(p - 1)c pif \mu \leqslant | \xi | (1 + c) \mu p - 1| \xi | - \mu p(p - 1) pif | \xi | (1 + c)< \mu < | \xi | 1 p| \xi | pif \mu \geqslant | \xi | . Observe that Fis nonstrictly convex in \xi , and it is not differentiable with respect to \mu . We also prove the convergence of the solutions of the discretized problem toward the solutions of the continuous one. Moreover, taking a regular sequence of triangulations in \Omega of diameter h > 0 and discretizing W1,p 0(\Omega ) and L\infty (\Omega ) by the usual P1and P0finite elements, respectively, we show that the difference between the minimum for the continuous and the discretized problem is of order h. In order to prove this result, we assume the existence of a solution (u, \theta ) for (1.3) such that uis in W1,\infty (\Omega ), \nabla u belongs to BV (\Omega )N, and \theta belongs to BV (\Omega ). Some smoothness results for problem (1.3) can be found in [6] and [14] for p= 2 and [7] for p\in (1,\infty ); we also refer the reader to [5] for the relaxed problem corresponding to take the minimum in (1.1) instead of the maximum one. These smoothness results imply that uis in W1,\infty (\Omega ), the flow \sigma =| \nabla u| p - 2\nabla u/(1 + c\theta )p - 1is in H1(\Omega )N, and the derivatives of \theta in the direction of \sigma are in L2(\Omega ). However, this is not enough to get \nabla uand \theta BV -functions. Nevertheless, this assumption seems to be satisfied in the numerical experiments. The paper is organized as follows. In section 2, we recall some known results for problem (1.3) which have been proved in [7] (see [6], [20], for p= 2). In section 3, we state the main results of the paper. Section 4 is devoted to prove the results in section 3. Finally, in section 5, we illustrate the results of the paper with some numerical simulations. They show that the alternating minimization method converges faster than the feasible direction method. 2. Previous results. As we mentioned in the introduction, our aim in the present paper is to numerically solve the optimal design problem (1.1). Since it has no solution in general, we work with the relaxed formulation (1.3), which, renaming f/\beta by fto simplify the notation, can be written as (2.1) min \biggl\{ \scrF (\theta , u) : \theta \in L\infty (\Omega ; [0,1]), u \in W1,p 0(\Omega ),\int \Omega \theta dx \leqslant \kappa \biggr\} , with (2.2) \scrF (\theta , u) = 1 p\int \Omega | \nabla u| p (1 + c \theta )p - 1dx - \langle f, u\rangle . Here \Omega is a bounded open set of \BbbR N,N\geqslant 2, p\in (1,\infty ), c > 0, \kappa \in (0,| \Omega | ), and fis a distribution in W - 1,p\prime (\Omega ). Since \scrF is convex in (\theta , u) and coercive in u,W1,p 0(\Omega ) is reflexive and L\infty (\Omega ; [0,1]) is bounded and then sequentially compact for the weak-\ast topology in L\infty (\Omega ), the existence of solution is straightforward. However, \scrF is not strictly convex, and therefore the uniqueness is not clear. Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 3080 J. CASADO-D\' IAZ, C. CONCA, AND D. V\' ASQUEZ-VARAS The relaxed formulation (2.1) has been obtained in [7]. In this paper, we have also obtained some optimality conditions and some equivalent formulations. As a consequence, we got some uniqueness and smoothness results (see [6], [12], [14], [20] for related results in the case p= 2). Thanks to the convexity of \scrF , Kuhn--Tucker's theorem easily provides the following system of optimality conditions ([7]). Proposition 2.1. A pair (\^ \theta , \^u)is a solution of (2.1) if and only if there exists \^\mu \geqslant 0such that the following hold: If \^\mu = 0, then (2.3) \^ \theta = 1 a.e. in \{ x\in \Omega : \nabla \^u(x)\not = 0\bigr\} ,\bigm| \bigm| \{ x\in \Omega : \nabla \^u(x)\not = 0\bigr\} \bigm| \bigm| \leqslant \kappa , (2.4) \left\{ - div \biggl( | \nabla \^u| p - 2 (1 + c)p - 1\nabla \^u\biggr) =fin \Omega , \^u= 0 on \partial \Omega . If \^\mu > 0, then (2.5) \^ \theta = max \biggl\{ 0,min \biggl\{ 1,1 c\biggl( | \nabla \^u| \^\mu - 1\biggr) \biggr\} \biggr\} ,\int \Omega \^ \theta dx =\kappa , (2.6) \left\{ - div \Biggl( | \nabla \^u| p - 2 (1 + c\^ \theta )p - 1\nabla \^u\Biggr) =fin \Omega , \^u= 0 on \partial \Omega . Remark 2.1. The expression of \^ \theta in Proposition 2.1 is obtained by solving (see [7], [20]) (2.7) min \biggl\{ \int \Omega | \nabla \^u| p (1 + c\theta )p - 1dx :\theta \in L\infty (\Omega ; [0,1]),\int \Omega \theta dx \leqslant \kappa \biggr\} . The constant \^\mu \geqslant 0 is a Lagrange multiplier corresponding to the constraint \int \Omega \theta dx \leqslant \kappa . We observe that for an arbitrary function \^u\in W1,p 0(\Omega ) (not necessarily a solution for (2.1)), the solutions of (2.7) can be explicitly obtained using Kuhn--Tucker's theorem, which shows that \^ \theta is a solution if and only if there exists \^\mu \geqslant 0 such that \^\mu \biggl( \int \Omega \^ \theta dx - \kappa \biggr) = 0 and \^ \theta is a solution of (2.8) min \biggl\{ \int \Omega | \nabla \^u| p (1 + c\theta )p - 1dx + \^\mu \int \Omega \^ \theta dx :\theta \in L\infty (\Omega ; [0,1])\biggr\} . This provides the following rule to solve (2.7). If \^uis such that | \{ \nabla \^u\not = 0\} | \leqslant \kappa , then \^ \theta is any function in L\infty (\Omega ; [0,1]) satisfying \^ \theta = 1 a.e. in \{ x\in \Omega : \nabla \^u(x)\not = 0\} ,\int \Omega \^ \theta dx \leqslant \kappa . In the other case, denoting for \mu > 0 \theta \mu := max \biggl\{ 0,min \biggl\{ 1,1 c\biggl( | \nabla \^u| \mu - 1\biggr) \biggr\} \biggr\} Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. NUMERICAL MAXIMIZATION OF THE p-LAPLACIAN ENERGY 3081 and defining G: (0,\infty )\rightarrow \bigl[ 0,| \Omega | \bigr] by G(\mu ) = \int \Omega \theta \mu dx \forall \mu \in (0,\infty ), we have that the set of solutions of (2.7) is given by (2.9) \bigl\{ \theta \mu \in L\infty (\Omega ; [0,1]) : \mu > 0, G(\mu ) = \kappa \bigr\} . Remark that the equation G(\mu ) = \kappa has a solution (not unique in general) due to G decreasing, continuous, and lim \mu \rightarrow 0G(\mu ) = \bigm| \bigm| \{ x\in \Omega : \nabla \^u(x)\not = 0\} \bigm| \bigm| ,lim \mu \rightarrow \infty G(\mu )=0. Numerically, the equation G(\mu ) = \kappa can be easily solved using, for example, a dichotomy method. In [7] (see [20] for p= 2), it has also been proved that introducing the flow \sigma =| \nabla u| p - 2 (1 + c\theta )p - 1\nabla u, we have that (2.1) is equivalent to the min-max problem (2.10) min - div \sigma =f \sigma \in Lp\prime (\Omega )N max \theta \in L\infty (\Omega :[0,1]) \int \Omega \theta dx\leqslant \kappa \int \Omega (1+c \theta )| \sigma | p\prime dx = max \theta \in L\infty (\Omega :[0,1]) \int \Omega \theta dx\leqslant \kappa min - div \sigma =f \sigma \in Lp\prime (\Omega )N\int \Omega (1+c \theta )| \sigma | p\prime dx. Taking into account that the functional \sigma \in L2(\Omega )N\rightarrow max \theta \in L\infty (\Omega :[0,1]) \int \Omega \theta dx\leqslant \kappa \int \Omega (1 + c \theta )| \sigma | p\prime dx is strictly convex, we get the uniqueness of the optimal flow. Moreover, using (1.4), we get the following smoothness results for the solutions of (2.1). Theorem 2.1. For every solution (\^ \theta , \^u)of (2.1), the flow \^\sigma defined by (2.11) \^\sigma =| \nabla \^u| p - 2 (1 + c\^ \theta )p - 1\nabla \^u is uniquely defined. If fbelongs to W1,1(\Omega ) \cap Lq(\Omega ),q > N, and \Omega is a C1,1domain, then \^\sigma belongs to H1(\Omega )N\cap L\infty (\Omega )N. Moreover, there exists C > 0, which only depends on N,p,c, and \Omega , such that (2.12) \| \^\sigma \| H1(\Omega )N\cap L\infty (\Omega )N\leqslant C\Bigl( \| f\| W1,1(\Omega )\cap Lq(\Omega ) + \^\mu \Bigr) , with \^\mu given by Proposition 2.1. The function \^ \theta satisfies (2.13) \^ \theta (x) = \Biggl\{ 1if | \^\sigma | >\^\mu , 0if | \^\sigma | <\^\mu , and decomposing \^\sigma = (\^\sigma 1,...,\^\sigma N), we have (2.14) \partial i\^ \theta \^\sigma j - \partial j\^ \theta \^\sigma i= (1 + c\^ \theta )(\partial j\^\sigma i - \partial i\^\sigma j)\chi \{ | \^\sigma | =\^\mu \} \in L2(\Omega ),1\leqslant i, j \leqslant N. Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 3082 J. CASADO-D\' IAZ, C. CONCA, AND D. V\' ASQUEZ-VARAS Remark 2.2. Theorem 2.1 has been proved in [7], where some other regularity results depending on the smoothness of fhave been obtained. The case p= 2 has been first proved in [6]. Observe that \^\sigma in L\infty (\Omega )Nimplies that \^ubelongs to W1,\infty (\Omega ). This was previously shown in [14] for p= 2. 3. Algorithms and main results. In this section, we present two variants of a descent algorithm to numerically solve problem (2.1). We also show the convergence of both algorithms. A first attempt to construct an algorithm is to use an alternate method consisting in minimizing in u, then in \theta , and so on. That is, assuming an approximation (ui, \theta i) of a solution of (2.1), we compute ui+1 as a solution of (3.1) min v\in W1,p 0(\Omega ) \biggl\{ 1 p\int \Omega | \nabla v| p (1 + c\theta i)p - 1 - \langle f, v\rangle \biggr\} and then \theta i+1 as a solution of (3.2) min \theta \in L\infty (\Omega ;[0,1]) \int \Omega \theta dx\leqslant \kappa \biggl\{ 1 p\int \Omega | \nabla ui+1| p (1 + c\theta )p - 1dx\biggr\} . This method works well if p= 2, but for p\not = 2, problem (3.1) is a p-Laplacian problem, which is very expensive to solve from the computational point of view due to the nonlinearity of the corresponding Euler--Lagrange equation. Instead of using the above alternate method, we can also try to use a gradient method, i.e., an iterative method where the iterations are defined through ui+1 = ui+tivi+1,\theta i+1 =\theta i+si(\vargamma i+1 - \theta i) for some ti, si\in (0,1), with (vi+1, \vargamma i+1) a solution of (3.3) \left\{ min \| v\| W1,p 0(\Omega )\leqslant 1\biggl\{ 1 p\int \Omega | \nabla ui| p - 2 (1 + c\theta i)p - 1\nabla ui\cdot \nabla v dx - \langle f, v\rangle \biggr\} , max \vargamma \in L\infty (\Omega ;[0,1]) \int \Omega \vargamma dx\leqslant \kappa \int \Omega | \nabla ui| p (1 + c\theta i)p\vargamma dx, but the minimization in valso implies the resolution of a p-Laplacian problem. To avoid this difficulty, we can replace the constraint \| v\| W1,p 0(\Omega ) \leqslant 1 by \| v\| H1 0(\Omega ) \leqslant 1. This is a feasible direction method. In each iteration, we look for the direction of maximum descent of \scrF in the convex set: \biggl\{ (v, \vargamma )\in H1 0(\Omega ) \times L\infty (\Omega ; [0,1]) : \| v\| H1 0(\Omega ) \leqslant 1,\int \Omega \vargamma dx \leqslant \kappa \biggr\} . The maximum direction with respect to \vargamma is simple to calculate. Namely, reasoning as in Remark 2.1, we have the following. If | \{ | \nabla ui| >0\} | \leqslant \kappa , then \vargamma iis any function in L\infty (\Omega ; [0,1]) such that (3.4) \chi \{ | \nabla ui| >0\} \leqslant \vargamma i,\int \Omega \vargamma idx \leqslant \kappa . In another case, we introduce H: (0,\infty )\rightarrow [0,| \Omega | ] by H(\mu ) = \bigm| \bigm| \{ x\in \Omega : | \nabla ui(x)| >(1 + c\theta i)\mu \bigr\} \forall \mu \geqslant 0. Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. NUMERICAL MAXIMIZATION OF THE p-LAPLACIAN ENERGY 3083 Then His a decreasing function, continuous on the right and satisfying lim \mu \rightarrow 0+H(\mu ) = \bigm| \bigm| \{ x\in \Omega : | \nabla ui(x)| >0\} \bigm| \bigm| ,lim \mu \rightarrow \infty H(\mu )=0. This ensures the existence of \mu i>0 (not unique in general) such that H(\mu i)\leqslant \kappa \leqslant lim \mu \rightarrow \mu - i H(\mu ), which can be easily numerically obtained by a dichotomy rule. For such \mu i, the maximum direction in \theta in (3.3), \vargamma i, is given by any function in L\infty (\Omega ; [0,1]) such that (3.5) \chi \{ | \nabla ui| >(1+c\theta i)\mu i\} \leqslant \vargamma i,\int \Omega \vargamma idx =\kappa . A similar result holds if we use a finite-dimensional approximation consisting in choosing \theta taking constant values in the elements of a given mesh. On the other hand, the maximum descent direction with respect to vis unique, and it is the solution of a linear equation. However, we observe that for p > 2, the sequence of functions \{ ui\} generated by the method is not in W1,p 0(\Omega ). Thus, the algorithm has only a sense of using a finite-dimensional space instead of L\infty (\Omega ) \times W1,p 0(\Omega ). In such a case, all the norms are equivalent. However it would be necessary to prove the convergence of the solutions of the discretized problem to the continuous one. With these considerations, we are going to be interested in the problem (3.6) min \biggl\{ \scrF (\theta , u) : \theta \in \Theta , u \in V, \int \Omega \theta dx \leqslant \kappa , \theta \in [0,1] a.e. in \Omega \biggr\} , with \Theta and Vfinite-dimensional subspaces of L\infty (\Omega ) and H1 0(\Omega ) \cap W1,p 0(\Omega ), respectively. As in the continuous problem, it is not clear that (3.6) has a unique solution, but for every solution (\theta \ast , u\ast ), the flow (3.7) \sigma \ast =| \nabla u\ast | p - 2 (1 + c\theta \ast )p - 1\nabla u\ast is unique because it is a solution of (see (2.10)) min \left\{ max \theta \in L\infty (\Omega :[0,1]) \int \Omega \theta dx\leqslant \kappa \int \Omega (1 + c \theta )| \sigma | p\prime dx :\sigma \in Lp\prime (\Omega )N,\int \Omega \sigma \cdot \nabla v dx =\langle f, v\rangle ,\forall v\in V\right\} , where the function to minimize is strictly convex. As an example of practical interest, we can consider a regular triangular mesh \scrT h of \Omega with maximum diameter h > 0 and the Lagrange finite element spaces \Theta h=\Bigl\{ v=\sum \tau \in \scrT h \alpha \tau \scrX \tau :\alpha \tau \in \BbbR \forall \tau \in \scrT h\Bigr\} ,(3.8) Vh=\bigl\{ v\in C0 0(\Omega ) : v| \tau \in \BbbP 1(\tau )\forall \tau \in \scrT h\bigr\} ,(3.9) with \BbbP 1(\tau ) the space of affine functions in \tau . Since the minimization of \scrF in \theta for ufixed is simple to carry out in practice (see (2.1) for the infinite-dimensional case; the finite-dimensional one is analogous), we Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 3084 J. CASADO-D\' IAZ, C. CONCA, AND D. V\' ASQUEZ-VARAS can also consider a variant of the previous algorithm consisting in directly computing the minimum in \theta in each iteration. With these considerations, we present the following two algorithms. Algorithm 1. Initialization: i= 1, \theta 0\in \Theta , u0\in V,a, b \in (0,1). 1: Set via solution of (3.10) \int \Omega \nabla vi\cdot \nabla \phi dx =\langle f, \phi \rangle - \int \Omega | \nabla ui| p - 2 (1 + c\theta i)p - 1\nabla ui\cdot \nabla \phi dx \forall \phi \in V. 2: Choose the step length by ti=bj(Armijo's rule), with jthe smallest nonnegative integer such that (3.11) \scrF (\theta i, ui+tivi)\leqslant \scrF (\theta i, ui) - ati\int \Omega | \nabla vi| 2dx, and set ui+1 =ui+tivi. 3: Set \vargamma ia solution of (3.12) max \biggl\{ \int \Omega | \nabla ui+1| p (1 + c\theta i)p\vargamma dx :\vargamma \in \Theta ,0\leqslant \vargamma \leqslant 1 a.e. in \Omega ,\int \Omega \vargamma dx \leqslant \kappa \biggr\} . 4: Choose si=bk, with kthe smallest nonnegative integer such that (3.13) \scrF (\theta i+si(\vargamma i - \theta i), ui+1)\leqslant \scrF (\theta i, ui+1) - asi c(p - 1) p\int \Omega | \nabla ui+1| p (1 + c\theta i)p(\vargamma i - \theta i)dx, and set \theta i+1 =\theta i+si(\vargamma i - \theta i). Algorithm 2. Initialization: i= 0, u0\in V,a, b \in (0,1). 1: Set vi\in Vthe solution of (3.10). 2: Choose the step length by ti=bjwith jthe smallest nonnegative integer such that (3.11) is satisfied, and set ui+1 =ui+tivi. 3: Set \theta i+1 a solution of (3.14) min \biggl\{ \int \Omega | \nabla ui| p (1 + c\vargamma )p - 1dx :\vargamma \in \Theta ,0\leqslant \vargamma \leqslant 1 a.e. in \Omega ,\int \Omega \vargamma dx \leqslant \kappa \biggr\} . Remark 3.1. Since by definition (3.10) of viwe have lim t\rightarrow 0 \scrF (\theta i, ui+tvi) - \scrF (\theta i, ui) t= - \int \Omega | \nabla vi| 2dx and a < 1, we get that \scrF (\theta i, ui+tvi) - \scrF (\theta i, ui)\leqslant - at \int \Omega | \nabla vi| 2dx, for 0 < t small enough. This proves the existence of tisatisfying (3.11). A similar argument shows the existence of siin (3.13). Remark 3.2. If p= 2, then ti= 1 for both algorithms and every i\geqslant 0. The second method agrees in this case with the one given in [1] Theorem 5.1.5, and [22]. Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. NUMERICAL MAXIMIZATION OF THE p-LAPLACIAN ENERGY 3085 Our main result is given by theorem 3.1 below which provides the convergence for both algorithms. Before stating it, we need the following definition. Definition 3.3. For p > 1, we define \gamma p>0by (3.15) \left\{ \| v\| H1 0(\Omega ) \leqslant \gamma p\| v\| W1,p 0(\Omega ) if 1<p<2 \| v\| W1.p 0(\Omega ) \leqslant \gamma p\| v\| H1 0(\Omega ) if p\geqslant 2\forall v\in V. Remark 3.4. Clearly \gamma 2= 1, while for p\not = 2 and Vreplaced by a sequence of finite-dimensional spaces Vhsuch that lim h\rightarrow 0min v\in Vh \| v - vh\| W1,p 0(\Omega ) = 0 \forall v\in W1.p 0(\Omega ), we have that \gamma p=\gamma p,h tends to infinity when hgoes to zero. For example, in the case where the spaces Vhare given by (3.9), with \scrT ha sequence of regular meshes of diameter h, we have (3.16) \gamma p,h =O\Biggl( 1 hN| 1 2 - 1 p| \Biggr) . Theorem 3.1. Let \Omega \subset \BbbR Nbe a bounded open set, p\in (1,\infty ),f\in W - 1,p\prime (\Omega ), and \Theta ,Vfinite-dimensional subspaces of L\infty (\Omega ) and W1,p 0(\Omega ), respectively. Taking (\theta i, ui)\in \Theta \times V, the sequence defined by Algorithm 1or 2, and denoting by \scrF \ast the minimum value of (3.6) and by eithe sequence of errors (3.17) ei=\scrF (\theta i, ui) - \scrF \ast \geqslant 0, i \geqslant 0, we have that eiis a decreasing sequence and that there exists C > 0depending on a, b, u0, \theta 0, f, c, N, and psuch that (3.18) ei\leqslant \Biggl\{ C\gamma p pi - 1 p - 1if 1<p<2 C\gamma 4 pi - 1if p\geqslant 2 \forall i\geqslant 1. Moreover, the sequence (3.19) \sigma i=| \nabla ui| p - 2 (1 + c\theta i)p - 1\nabla ui converges strongly to \sigma \ast defined by (3.7) in Lp\prime (\Omega )N. Namely, there exists C > 0as above such that (3.20) \int \Omega \bigl( | \sigma \ast | +| \sigma i| )p - 2| \sigma \ast - \sigma | 2dx \leqslant \left\{ C\bigl( 1 + \gamma p\| u0\| W1,p 0(\Omega )\bigr) (ei - ei+1)1 p\prime if 1<p<2, C\bigl( | 1 + \gamma 2 p\| u0\| p - 1 W1,p 0(\Omega )\bigr) (ei - ei+1)1 2if p\geqslant 2. Remark 3.5. In the continuous case V=W1,p 0(\Omega ) and 1 < p < 2, the classical regularity results for the Poisson equation show that the solution viof (3.10) satisfies the estimate (3.21) \| vi\| W1,p\prime 0(\Omega ) \leqslant C\bigl( \| ui\| W1,p 0(\Omega ) +\| f\| W - 1,p\prime (\Omega )\bigr) , with C > 0 depending only on pand \Omega . Thanks to this result, we can deduce that in this case (3.18) holds true with \gamma preplaced by one. Similar results to (3.21) also hold for special choices of spaces V; see, e.g., [4, Theorem 8.5.3]. With these choices, we can eliminate the dependence in \gamma pof estimate (3.18) for 1 <p<2. Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 3092 J. CASADO-D\' IAZ, C. CONCA, AND D. V\' ASQUEZ-VARAS with \rho positive constants which only depend on p. Similarly, using \vargamma ias a solution of (3.12), we have (4.20) \int \Omega (1 + c\vargamma i)| \sigma i| p\prime dx \geqslant \int \Omega (1 + c\theta \ast )| \sigma \ast | p\prime dx +p\prime \int \Omega (1 + c\theta \ast )| \sigma \ast | p\prime - 2\sigma \ast \cdot (\sigma i - \sigma \ast )dx +\rho \int \Omega \bigl( | \sigma \ast | +| \sigma i| )p - 2| \sigma \ast - \sigma i| 2dx. From (4.19) and (4.20), we deduce (4.21) 0\geqslant p\prime \int \Omega \Bigl( (1 + c\theta i)| \sigma i| p\prime - 2\sigma i - (1 + c\theta \ast )| \sigma \ast | p\prime - 2\sigma \ast \Bigr) \cdot (\sigma \ast - \sigma i)dx + 2\rho \int \Omega \bigl( | \sigma \ast | +| \sigma i| )p - 2| \sigma \ast - \sigma i| 2dx +c\int \Omega (\theta i - \vargamma i)| \sigma i| p\prime dx. Now we use that \Bigl( (1 + c\theta i)| \sigma i| p\prime - 2\sigma i - (1 + c\theta \ast )| \sigma \ast | p\prime - 2\sigma \ast \Bigr) \cdot (\sigma \ast - \sigma i) =\Bigl( | \nabla u\ast | p - 2 (1 + c\theta \ast )p - 1\nabla u\ast - | \nabla ui| p - 2 (1 + c\theta i)p - 1\nabla ui\Bigr) \cdot \nabla (ui - u\ast ), which, taking into account (3.10) and that (\theta \ast , \sigma \ast ) satisfies the discrete version of (2.7), prove (4.22) \int \Omega \Bigl( (1 + c\theta i)| \sigma i| p\prime - 2\sigma i - (1 + c\theta \ast )| \sigma \ast | p\prime - 2\sigma \ast \Bigr) \cdot (\sigma \ast - \sigma i)dx =\int \Omega \nabla (ui - u\ast )\cdot \nabla vidx. Replacing this equality in (4.21) and recalling ui+1 =ui+tivi, with 0 \leqslant ti\leqslant 1, we get (4.23) 2\rho \int \Omega \bigl( | \sigma \ast | +| \sigma i| )p\prime - 2| \sigma \ast - \sigma i| 2dx \leqslant c\int \Omega | \nabla ui+1| p (1 + c\theta i)p(\vargamma i - \theta i)dx +C\int \Omega \bigl( | \nabla ui| +| \nabla u\ast | \bigr) p - 1| \nabla vi| dx +p\prime \int \Omega \bigl( | \nabla ui| +| \nabla u\ast | \bigr) | \nabla vi| dx with Cdepending only on pand c. By (4.13) and (4.14), we then conclude (3.20). Proof of Theorem 3.2.For (\^ \theta , \^u) the solution of (2.1) which satisfies (3.28), \^\sigma defined by (2.11), and h > 0, we introduce \^ \theta h\in \Theta h, \^\sigma h\in \Theta N h, and \^uh\in Vhby (4.24) \^ \theta h=1 | \tau | \int \tau \^ \theta dx, \^\sigma h=1 | \tau | \int \tau \^\sigma dx \forall \tau \in \scrT h, (4.25) \^uh(xi) = \^u(xi)\forall xivertex of \scrT h. Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. NUMERICAL MAXIMIZATION OF THE p-LAPLACIAN ENERGY 3093 Thanks to (3.28) and the regularity of \scrT h, there exists C > 0 such that (4.26) h\| \^uh\| W1,\infty (\Omega ) +\| \^uh - \^u\| L\infty (\Omega ) +\| \^uh - \^u\| W1,1 0(\Omega ) +\| \^ \theta h - \^ \theta \| L1(\Omega ) \leqslant Ch. The definition of \scrF , the mean value theorem, and these estimates imply \bigm| \bigm| \bigm| \scrF (\^ \theta , \^u) - \scrF (\^ \theta h,\^uh)\bigm| \bigm| \bigm| \leqslant C\Bigl( \| \nabla \^u\| L\infty (\Omega )N+\| \nabla \^uh\| L\infty (\Omega )N\Bigr) p - 1\| \nabla (u - uh)\| L1(\Omega )N +\| \nabla \^u\| p L\infty (\Omega )N\| \^ \theta - \^ \theta h\| L1(\Omega ) +\| f\| L1(\Omega )\| \^u - \^uh\| L\infty (\Omega ) \leqslant Ch. Then, since the definitions of \^ \scrF and \scrF \ast himply \scrF (\^ \theta , \^u) = \^ \scrF \leqslant \scrF \ast h\leqslant \scrF (\^ \theta h,\^uh), we conclude (3.29). On the other hand, we consider (\theta \ast h, u\ast h) a solution of (3.6) with \Theta and Vreplaced by \Theta hand Vh. We define \^\sigma by (2.11) and \sigma \ast hby (4.27) \sigma \ast h=| \nabla u\ast h| p - 2 (1 + c\theta \ast h)p - 1\nabla u\ast h, and we recall that, thanks to (2.10), we have \int \Omega (1 + c\^ \theta )| \^\sigma | p\prime dx = max \biggl\{ \int \Omega (1 + c\theta )| \^\sigma | p\prime dx :\theta \in L\infty (\Omega ; [0,1]),\int \Omega \theta dx \leqslant \kappa \biggr\} . Thus, we deduce (4.28) \int \Omega (1 + c\^ \theta )| \^\sigma | p\prime dx \geqslant \int \Omega (1 + c\theta \ast h)| \^\sigma | p\prime dx \geqslant \int \Omega (1 + c\theta \ast h)| \sigma \ast h| p\prime dx +p\prime \int \Omega (1 + c\theta \ast h)| \^\sigma \ast h| p\prime - 2\^\sigma \ast h\cdot (\^\sigma - \^\sigma \ast h)dx +\rho \int \Omega (| \^\sigma | + \^\sigma \ast h| )p - 2| \^\sigma - \^\sigma \ast h| 2dx for some \rho > 0, which only depends on p. Using the definitions of \^\sigma and \^\sigma \ast and that (\^ \theta , \^u), (\theta \ast , u\ast h) are solutions of (2.1) and (3.6), we have \int \Omega (1 + c\^ \theta )| \^\sigma | p\prime dx =\int \Omega | \nabla \^u| p (1 + c\^ \theta )p - 1dx = - p\prime \^ \scrF , \int \Omega (1 + c\theta \ast h)| \sigma \ast h| p\prime dx =\int \Omega | \nabla u\ast h| p (1 + c\theta \ast h)p - 1dx = - p\prime \scrF \ast h, p\prime \int \Omega (1 + c\theta \ast h)| \^\sigma \ast h| p\prime - 2\^\sigma \ast h\cdot (\^\sigma - \^\sigma \ast h)dx =p\prime \int \Omega (\^\sigma - \^\sigma \ast h)\cdot \nabla u\ast hdx = 0. Replacing these equalities in (4.28) and taking into account (3.29), we get (3.30). Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 3094 J. CASADO-D\' IAZ, C. CONCA, AND D. V\' ASQUEZ-VARAS 5. Numerical experiments. In this section, we present some simulations for the numerical resolution of (2.1) using the two algorithms presented in section 3. The implementation has been carried out in Python using the finite element solver FeniCs [2]. In our numerical experiments, we have taken N= 2, \Omega the unit disc, c= 1, f = 1, and \kappa = 1. In this case, the solution of (2.1) is explicitly given by \^ \theta (x) = \Biggl\{ 1 if | x| < \pi - 1 2, 0 if | x| > \pi - 1 2, \^u(x) = \left\{ 1 21 p - 1p\prime \bigl( 1 + \pi - p\prime 2 - 2| x| p\prime \bigr) if | x| < \pi - 1 2, 1 21 p - 1p\prime \bigl( 1 - | x| p\prime \bigr) if | x| > \pi - 1 2, and thus \^ \scrF =\scrF (\^ \theta , \^u) = - \pi p\prime (2 + p\prime )2 1 p - 1\Bigl( 1 - 1 2\pi 1+ p\prime 2\Bigr) . We solve the problem for meshes of different diameter hand p= 1.2,2,100. The stop criterion for the first algorithm is (5.1) \int \Omega | \nabla vi| 2dx +c(p - 1) p\int \Omega | \nabla ui+1| p (1 + \theta i)p(\vargamma i - \theta i)dx \leqslant 10 - 7or i\geqslant 2000, while for the second one, it is given by (5.2) \int \Omega | \nabla vi| 2dx +c(p - 1) p\int \Omega | \nabla ui+1| p (1 + \theta i+1)p(\theta i+1 - \theta i)dx \leqslant 10 - 7or i\geqslant 2000. Observe that in both cases, replacing 10 - 7by 0 would mean that (\theta i, ui) satisfies the optimality conditions for (2.1) and then, by the convexity of \scrF , that (\theta i, ui) is a solution for (3.6). Depending on p,h, and the choice of the algorithm, we present in Figure 1 the convergence history of the objective function, the Lagrange multiplier \mu , and the stop criterion (\| DF\| denotes the left-hand sides in (5.1) and (5.2), respectively). Observe that for p= 1.2 and p= 2, the rate of convergence for both algorithms is similar. However, for p= 100, Algorithm 2 converges faster than Algorithm 1. Although our estimates depend on h, we do not observe this dependence in the numerical experiments for p= 2. This is because \gamma 2= 1, and therefore, according to Remark 3.2, the step length is constant, and all the bounds in Theorem 3.1 do not depend on the mesh size. In Figure 2, we represent the solutions (\theta i, ui) depending on pbut only for the finest mesh. Observe that the solutions obtained are very similar for both algorithms. In Figure 3, we show the time spent in the resolution of the numerical experiments. We observe that the iterations are calculated faster for Algorithm 2 than for Algorithm 1. When the diameter of the mesh decreases, the time increases for both algorithms in the same way. Moreover, for plarge, Algorithm 2 needs fewer iterations than Algorithm 1, while for psmall, both algorithms use more or less the same number of iterations. In Figure 4, we present the value of the objective function for the final interaction for each mesh size and Algorithm 1 and 2. For p= 1.2 and p= 2, both algorithms have the same behaviours. On the other hand, when p= 100 Algorithm 2 achieves a lower value of the objective function in the final iteration for all of the mesh sizes. This agrees with the fact that Algorithm 2 is faster than Algorithm 1, as illustrated in Figure 2. Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. NUMERICAL MAXIMIZATION OF THE p-LAPLACIAN ENERGY 3095 Fig. 1.Convergence history for each pand mesh. Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 3096 J. CASADO-D\' IAZ, C. CONCA, AND D. V\' ASQUEZ-VARAS Fig. 2.Solutions for the finest mesh. Fig. 3.Mean CPU time by iteration in seconds for each p. Fig. 4.Convergence rate of minimum value as function of the mesh size for each p. Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. NUMERICAL MAXIMIZATION OF THE p-LAPLACIAN ENERGY 3097 REFERENCES [1] G. Allaire,Shape Optimization by the Homogenization Method, Appl. Math. Sci. 146, Springer-Verlag, Berlin, 2002. [2] M. S. Aln{\ae}s, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells,The FeniCs Project Version 1.5, Arch. Numer. Softw., 3 (2015), https://doi.org/10.11588/ans.2015.100.20553. [3] L. Blank and C. Rupprecht,An extension of the projected gradient method to a Banach space setting with application in structural topology optimization, SIAM J. Control Optim., 55 (2017), pp. 1481--1499. [4] S. C. Brenner and L. R. Scott,The Mathematical Theory of Finite Element Methods, 3rd ed., Text Appl. Math. 15, Springer-Verlag, Berlin, 2008. [5] J. Casado-D\' {\i}az,Some smoothness results for the optimal design of a two-composite material which minimizes the energy, Calc. Var. Partial Differential Equations, 53 (2015), pp. 649-- 673. [6] J. Casado-D\' {\i}az,Smoothness properties for the optimal mixture of two isotropic materials: The compliance and eigenvalue problems, SIAM J. Control Optim., 53 (2015), pp. 2319--2349. [7] J. Casado-D\' {\i}az, C. Conca, and D. V\' asquez-Varas,The maximization of the p-Laplacian energy for a two-phase material, SIAM J. Control Optim., 59 (2021), pp. 1497--1519. [8] C. Conca, A. Laurain, and R. Mahadevan,Minimization of the ground state for two phase conductors in low contrast regime, SIAM J. Appl. Math., 72 (2012), pp. 1238--1259. [9] B. Kawohl,On a family of torsional creep problems, J. Reine Angew. Math., 410 (1990), pp. 1--22. [10] R. Glowinski and A. Marrocco,Sur l'approximation, par \'el\'ements finis d'ordre un, et la r\'esolution, par p\'enalisation-dualit\'e d'une classe de probl\`emes de Dirichlet non lin\'eaires, Rev. Fran\c caise Automat. Informat. Rech. Op\'er. S\'er. Rouge Anal. Num\'er., 9 (1975), pp. 41-- 76. [11] S. Gonz\' alez-Andrade,A preconditioned descent algorithm for variational inequalities of the second kind involving the p-Laplacian operator, Comput. Optim. Appl., 66 (2017), pp. 123-- 162. [12] J. Goodman, R. V. Kohn, and L. Reyna,Numerical study of a relaxed variational problem from optimal design, Comput. Methods Appl. Mech. Engrg., 57 (1986), pp. 107--127. [13] Y. Q. Huang, R. Li, and W. Liu,Preconditioned descent algorithms for p-Laplacian, J. Sci. Comput., 32 (2007), pp. 343--371. [14] B. Kawohl, J. Stara, and G. Wittum,Analysis and numerical studies of a problem of shape design, Arch. Ration. Mech. Anal., 114 (1991), pp. 343--363. [15] A. Laurain,Global minimizer of the ground state for two phase conductors in low contrast regime, ESAIM Control Optim. Calc. Var., 20 (2014), pp. 362--388. [16] A. Logg,Efficient representation of computational meshes, Int. J. Comput. Sci. Eng., 4 (2009), pp. 283--295. [17] F. Murat,Un contre-exemple pour le probl\`eme du contr\^ole dans les coefficients, C. R. Acad. Sci. Paris A, 273 (1971), pp. 708--711. [18] F. Murat, Th\'eor\`emes de non existence pour des probl\`emes de contr\^ole dans les coefficients, C. R. Acad. Sci. Paris A, 274 (1972), pp. 395--398. [19] F. Murat,H-convergence, S\'eminaire d'Analyse Fonctionnelle et Num\'erique, 1977--78, Universit\'e d'Alger, multicopied, 34 pp.; English translation, F. Murat and L. Tartar, H- convergence, Topics in the Mathematical Modelling of Composite Materials, L. Cherkaev and R. V. Kohn, eds., Birkha\"user, Boston, 1998, pp. 21--43. [20] F. Murat and L. Tartar,Calcul des variations et homog\'en\'eisation, Les m\'ethodes de l'homog\'en\'eisation: Th\'eorie et applications en physique, Eyrolles, Paris, 1985, pp. 319-- 369; English translation, F. Murat and L. Tartar, Calculus of Variations and Homogenization, Topics in the Mathematical Modelling of Composite Materials, L. Cherkaev and R. V. Kohn, eds., Birkha\"user, Boston, 1998. [21] L. Tartar,The General Theory of Homogenization: A Personalized Introduction, Springer- Verlag, Berlin, 2009. [22] A. M. Toader,The convergence of an algorithm in numerical shape optimization, C. R. Acad. Sci. Paris S\'er. I Math., 323 (1997), pp. 195--198. Downloaded 12/09/22 to 150.214.182.233 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy