\documentclass[reqno]{amsart} \AtBeginDocument{{\noindent\small {\em Electronic Journal of Differential Equations}, Vol. 2004(2004), No. 75, pp. 1--26.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu (login: ftp)} \thanks{\copyright 2004 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2004/75\hfil Constructive Sobolev gradient preconditioning] {Constructive Sobolev gradient preconditioning for semilinear elliptic systems} \author[J\'anos Kar\'atson \hfil EJDE-2004/75\hfilneg] {J\'anos Kar\'atson} % in alphabetical order \address{Department of Applied Analysis\\ ELTE University, H-1518 Budapest, Hungary} \email{karatson@cs.elte.hu} \date{} \thanks{Submitted March 18, 2004. Published May 21, 2004.} \thanks{Partially supported by the Hungarian National Research Fund OTKA} \subjclass[2000]{35J65, 49M10} \keywords{Sobolev gradient, semilinear elliptic systems, numerical solution, \hfill\break\indent preconditioning} \begin{abstract} We present a Sobolev gradient type preconditioning for iterative methods used in solving second order semilinear elliptic systems; the $n$-tuple of independent Laplacians acts as a preconditioning operator in Sobolev spaces. The theoretical iteration is done at the continuous level, providing a linearization approach that reduces the original problem to a system of linear Poisson equations. The method obtained preserves linear convergence when a polynomial growth of the lower order reaction type terms is involved. For the proof of linear convergence for systems with mixed boundary conditions, we use suitable energy spaces. We use Sobolev embedding estimates in the construction of the exact algorithm. The numerical implementation has focus on a direct and elementary realization, for which a detailed discussion and some examples are given. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{remark}[theorem]{Remark} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{corollary}[theorem]{Corollary} \newtheorem{lemma}[theorem]{Lemma} \newtheorem{definition}[theorem]{Definition} \newtheorem{example}[theorem]{Example} \section{Introduction} The numerical solution of elliptic problems is a topic of basic importance in numerical mathematics. It has been a subject of extensive investigation in the past decades, thus having vast literature (cf. \cite{Ax-Lay, Hack, Kelley, O-R} and the references there). The most widespread way of finding numerical solutions is first discretizing the elliptic problem, then solving the arising system of algebraic equations by a solver which is generally some iterative method. In the case of nonlinear problems, most often Newton's method is used. However, when the work of compiling the Jacobians exceeds the advantage of quadratic convergence, one may prefer gradient type iterations including steepest descent or conjugate gradients (see e.g. \cite{Ax-G, Da}). An important example in this respect is the Sobolev gradient technique, which represents a general approach relying on descent methods and has provided various efficient numerical results \cite{NeuS, Neu1, Neu3}. In the context of gradient type iterations the crucial point is most often preconditioning. Namely, the condition number of the Jacobians of the discretized systems tends to infinity when discretization is refined, hence suitable nonlinear preconditioning technique has to be used to achieve a convenient condition number \cite{Ax, Ax-C}. The Sobolev gradient technique presents a general efficient preconditioning approach where the preconditioners are derived from the representation of the Sobolev inner product. The Sobolev gradient idea does in fact opposite to that which first discretizes the problem. Namely, the iteration may be theoretically defined in Sobolev spaces for the boundary-value problem (i.e. at the continuous level), reducing the nonlinear problem to auxiliary linear Poisson equations. Then discretization may be used for these auxiliary problems. This approach is based on the various infinite-dimensional generalizations of iterative methods, beginning with Kantorovich. For recent and earlier results see \cite{GGZ, K-A, Neu1, Neu3}. The author's investigations include the development of the preconditioning operator idea as shown in \cite{book}. Some recent numerical results are given in \cite{varprec, elpot} which are closely related to the Sobolev gradient idea. Namely, a suitable representation of the gradient yields a preconditioning elliptic operator; for Dirichlet problems the usual Sobolev inner product leads to the Laplacian as preconditioner. For systems one may define independent Laplacians as preconditioners, see \cite{JAA} for an earlier treatment for uniformly elliptic Dirichlet problems using the strong form of the operators. We note that the constructive representation of the Sobolev gradients with Laplacians in \cite{JAA, elpot} is due to a suitable regularity property. In this context the Sobolev gradient can be regarded as infinite dimensional preconditioning by the Laplacian. It yields that the speed of linear convergence is determined by the ellipticity properties of the original problem instead of the discretized equation, i.e., the ratio of convergence is explicitly obtained from the coefficients of the boundary-value problem. Therefore, it is independent of the numerical method used for the auxiliary problems. Another favourable property is the reduction of computational issues to those arising for the linear auxiliary problems. These advantages appear in the finite element methods (FEM) realization \cite{GFEM}. In \cite{GFEM, JAA, elpot} Dirichlet problems are considered for uniformly elliptic equations. The aim of this paper is to develop the above described realization of Sobolev gradients for semilinear elliptic systems, including the treatment of non-uniformity (caused by polynomial growth of the lower order reaction-type terms) such that linear convergence is preserved, and considering general mixed boundary conditions which need suitable energy spaces. The studied class of problems includes elliptic reaction-diffusion systems related to reactions of autocatalytic type. The paper first gives a general development of the method: after a brief Hilbert space background, the theoretical iteration is constructed in Sobolev space and convergence is proved. Linear convergence is obtained in the norm of the corresponding energy space, equivalent to the Sobolev norm. An excursion to Sobolev embeddings is enclosed, which is necessary for determining the required descent parameter (and the corresponding convergence quotient). Then numerical realization is considered with focus on direct elementary realization. A detailed discussion is devoted to the latter, giving a general extension of the ideas of \cite{nonloc, elpot}. Also numerical examples are presented. \section{General construction and convergence} % 2 \label{secgen} \subsection{The gradient method in Hilbert space} In this subsection the Hilbert space result of \cite{JAA} on non-differentiable operators is suitably modified for our purposes. First we quote the theorem on differentiable operators this result relies on. The classical theorem, mentioned already by Kantorovich \cite{K-A}, is given in the form needed for our setting, including suitable conditions and stepsize. \begin{theorem} \label{th21} Let $H$ be a real Hilbert space and $F :H\to H$ have the following properties: \begin{enumerate} \item[\rm (i)] $F$ is G\^ateaux differentiable; \item[\rm (ii)] for any $u,k,w,h\in H$ the mapping $s,t\mapsto F'(u+sk+tw)h$ is continuous from ${\bf R}^2$ to $H$; \item[\rm (iii)] for any $u\in H$ the operator $F'(u)$ is self-adjoint; \item[\rm (iv)] there are constants $M\ge m>0$ such that for all $u\in H$ \end{enumerate} $$m\|h\|^2 \le \langle F'(u)h,h\rangle \le M\|h\|^2 \quad (h\in H). $$ Then for any $b\in H$ the equation $F(u)=b$ has a unique solution $u^\ast \in H$, and for any $u^0\in H$ the sequence $$u^{k+1}= u^k -\frac{2}{M+m}(F(u^k)-b)\quad (k\in\mathbb{N}) $$ converges linearly to $ u^\ast$, namely, \begin{equation} \| u^k - u^\ast \| \le \frac{1}{m}\|F(u^0)-b \| \big( \frac{M-m}{M+m}\big)^k \quad (k\in\mathbb{N})\, . \label{B1} \end{equation} \end{theorem} A short proof of the theorem (cf. \cite{Ne}) is based on proving the estimate \begin{equation} \|F(u^k)-b \| \le \big( \frac{M-m}{M+m}\big)^k \|F(u^0)-b \| \quad (k\in\mathbb{N}) \label{Pe} \end{equation} (to which end one verifies that $J(u)\equiv u -\frac{2}{M+m}(F(u^k)-b)$ possesses contraction constant $\frac{M-m}{M+m}$). Then (\ref{Pe}) yields (\ref{B1}) since assumption (iv) implies \, $$ m\|u-v\|\le \|F(u)-F(v)\| \quad (u,v\in H). $$ \begin{proposition} \label{prop21} Under the assumptions of Theorem \ref{th21} we have $$ \|u^k-u^0\| < {1 \over m} \|F(u^0)-b \| \quad (k\in\mathbb{N}). $$ \end{proposition} \begin{proof} \begin{align*} \|u^k-u^0\| &\le \sum_{i=0}^{k-1}\|u_{i+1}-u_i\| \\ &= \frac{2}{M+m} \sum_{i=0}^{k-1} \|F(u^i)-b\| \\ &\le \frac{2}{M+m} \|F(u^0)-b\| \sum_{i=0}^{k-1} \big( \frac{M-m}{M+m}\big )^i \\ &\le \|F(u^0)-b\| \frac{2}{(M+m)(1-\frac{M-m}{M+m})} = {1 \over m} \|F(u^0)-b \|. \end{align*} \end{proof} Proposition \ref{prop21} allows localization of the ellipticity assumption (cf. \cite{GGZ}): \begin{corollary} \label{cor21} Let $u^0\in H$, $r_0= {1 \over m} \|F(u^0)-b \|$, $B(u^0,r_0)=\{u\in H:\, \|u-u^0\|\le r_0\}$. Then in Theorem \ref{th21} it suffices to assume that \begin{enumerate} \item[\rm (iv)'] there exist $M\ge m>0$ such that for all $u\in H$ $ \langle F'(u)h,h \rangle \ge m\|h\|^2 \quad (h\in H)$, and for all $u\in B(u^0,r_0)$ $\langle F'(u)h,h \rangle \le M\|h\|^2 \quad (h\in H)$, \end{enumerate} instead of assumption (iv), in order that the theorem holds. \end{corollary} \begin{proof} The lower bound in (iv)' ensures that $F$ is uniformly monotone, which yields existence and uniqueness as before. Owing to Proposition \ref{prop21} the upper bound $M$ is only exploited in $B(u^0,r_0)$ to produce the convergence result. \end{proof} Turning to non-differentiable operators, we now quote the corresponding result in \cite{JAA}. First the necessary notations are given. \begin{definition} \label{def21} \rm Let $B:D\to H$ be a strictly positive symmetric linear operator. Then the energy space of $B$, i.e. the completion of $D$ with respect to the scalar product $$ \langle x,y \rangle _B\equiv \langle Bx,y \rangle \quad (x,y\in D) $$ is denoted by $H_B$. The corresponding norm is denoted by $\|\cdot\|_B$. \end{definition} For any $r\in\mathbb{N}^+$ we denote by $H^r \equiv H\times H\times \dots \times H$ ($r$ times) the product space. The corresponding norm is denoted by $$ [[ u ]] \equiv \Big( \sum_{i=1}^{r} \| u_i \|^2 \Big)^{1/2} \quad (u\in H^r). $$ The obvious analogous notation is used for the products of $H_B$. \begin{theorem}[\cite{JAA}] \label{th22} Let $H$ be a real Hilbert space, $D\subset H$. Let $T_i:D^r\to H \quad (i=1,\dots ,r)$ be non-differentiable operators. We consider the system \begin{equation} T_i(u_1,\dots ,u_r)=g_i \quad (i=1,\dots ,r) \label{sy22} \end{equation} with $g=(g_1,\dots ,g_r)\in H^r$. Let $B:D\to H$ be a symmetric linear operator with lower bound $\lambda>0$. Assume that the following conditions hold: \begin{enumerate} \item[\rm (i)] $R(B) \supset R(T_i) \quad (i=1,\dots ,r)$; \item[\rm (ii)] for any i=1,\dots ,r the operators $B^{-1} T_i$ have G\^ateaux differentiable extensions $F_i:H_B^r\rightarrow H_B$ , respectively; \item[\rm (iii)] for any $u,k,w,h\in H_B^r$ the mappings $s,t \mapsto F_i'(u+sk+tw)h$ are continuous from ${\bf R}^2$ to $H_B$; \item[\rm (iv)] for any $u,h,k\in H_B^r$ $$ \sum_{i=1}^{r} \langle F_i'(u)h,k_i \rangle _B = \sum_{i=1}^{r}\langle h_i,F_i'(u)k \rangle _B ; $$ \item[\rm (v)] there are constants $M\ge m>0$ such that for all $u,h\in H_B^r$ $$ m\sum_{i=1}^{r}\|h_i\|^2_B \le \sum_{i=1}^{r} \langle F_i'(u)h, h_i \rangle_B \le M \sum_{i=1}^{r} \|h_i\|^2_B\, . $$ \end{enumerate} Let $g_i \in R(B) \quad (i=1,\dots ,r)$. Then \begin{enumerate} \item[\rm (1)] system (\ref{sy22}) has a unique weak solution $u^\ast=(u^\ast_1,\dots ,u^\ast_r) \in H_B^r$, i.e. which satisfies $$ \langle F_i(u^\ast ),v \rangle _B= \langle g_i,v \rangle \quad (v\in H_B,\, i=1,\dots ,r) ; $$ \item[\rm (2)] for any $u^0\in D^r$ the sequence $u^k=(u^k_1,\dots ,u^k_r)_{k\in\mathbb{N}}$, given by the coordinate sequences $$ u^{k+1}_i\equiv u^k_i - \frac{2}{M+m} B^{-1}(T_i(u^k)-g_i) \quad (i=1,\dots ,r; k\in\mathbb{N}), $$ converges linearly to $u^\ast $. Namely, $$ [[ u^k - u^\ast ]]_B \le \frac{1}{m\sqrt{\lambda}} [[ T(u^0)-g ]] \big( \frac{M-m}{M+m}\big) ^k \quad (k\in\mathbb{N})\, . $$ \end{enumerate} \end{theorem} The proof of this theorem in \cite{JAA} is done by applying Theorem \ref{th21} to the operator $F=(F_1,\dots ,F_r)$ and the right-side $b=\{ b_i\}_{i=1}^r =\{ B^{-1}g_i \}_{i=1}^r$ in the space $H_B^r$. This implies that the assumption can be weakened in the same way as in Corollary \ref{cor21}. That is, we have \begin{corollary} \label{cor22} Let $u^0\in D$, $b_i=B^{-1}g_i$ $(i=1,\dots ,r)$, $r_0= {1 \over m} [[F(u^0)-b ]]_B$, $B(u^0,r_0)=\{u\in H_B^r:\, [[u-u^0]]_B \le r_0 \}$. Then in Theorem 2.2 it suffices to assume that \begin{enumerate} \item[\rm (v)'] there exist $M\ge m>0$ such that for all $u\in H_B$ $ \sum_{i=1}^{r} \langle F_i'(u)h, h_i \rangle_B \ge m[[h]]_B^2 \quad (h\in H_B)$, and for all $u\in B(u^0,r_0)$ \ $ \sum_{i=1}^{r} \langle F_i'(u)h, h_i \rangle_B \le M[[h]]_B^2 \quad (h\in H_B)$, \end{enumerate} instead of assumption (v), in order that the theorem holds. \end{corollary} Finally we remark that the {\it conjugate gradient method (CGM)} in Hilbert space is formulated in \cite{Da} for differentiable operators under fairly similar conditions as for Corollary \ref{cor21}, and this result is extended to non-differentiable operators in \cite{cgmn} similarly to Corollary \ref{cor22}. Compared to the gradient method, the CGM improves the above convergence ratio to $(\sqrt{M}-\sqrt{m})(\sqrt{M}+\sqrt{m})$, on the other hand, the extra work is the similar construction of two simultaneous sequences together with the calculation of required inner products and numerical root finding for the stepsize. \subsection{The gradient method in Sobolev space} We consider the system of boundary value problems \begin{equation} \label{SY} \begin{gathered} T_i(u_1,\dots ,u_r)\equiv\, -\mathop{\rm div} (a_i(x)\nabla u_i) + f_i(x,u_1,\dots ,u_r )= g_i(x) \quad \mbox{in } \Omega \\ Qu_i \equiv (\alpha(x) u_i +\beta(x) \partial_\nu u_i) {}\big|_{\partial \Omega }=0 \end{gathered} \end{equation} $(i=1,\dots ,r)$ on a bounded domain $\Omega\subset\mathbb{R}^N$ with the following conditions: \begin{enumerate} \item[\rm (C1)] $\partial\Omega\in C^2$, $a_i\in C^1(\overline\Omega)$, $f_i\in C^1(\overline\Omega\times\mathbb{R}^r)$, $g_i\in L^2(\Omega)$. \item[\rm (C2)] $\alpha,\beta\in C^1(\partial\Omega)$, $\alpha,\beta\ge 0$, $\alpha^2 +\beta^2 >0$ almost everywhere on $\partial\Omega$. \item[\rm (C3)] There are constants $m,m'>0$ such that $00 \}. $$ \item[\rm (C4)] Let $2\le p\le {2N \over N-2}$ \, (if $N>2$), \, $2\le p$\, (if $N=2$). There exist constants $\kappa'\ge \kappa\ge 0$ and $\gamma\ge 0$ such that for any $(x,\xi)\in \overline\Omega\times\mathbb{R}^r$ the Jacobians $ \partial_\xi f (x,\xi) = \left\{ \partial_{\xi_k} f_j (x,\xi_1,\dots ,\xi_r) \right\}_{j,k=1}^r \in \mathbb{R}^{r\times r}$ are symmetric and their eigenvalues $\mu$ fulfil $$ \kappa \le \mu \le \kappa'+\gamma \sum_{j=1}^r |\xi_j|^{p-2}. $$ Moreover, in the case $\alpha\equiv 0$ we assume $\kappa >0$, otherwise $\kappa=0$. \end{enumerate} Let \begin{equation} D_Q\equiv \{ u\in H^2(\Omega): \, Qu\big|_{\partial\Omega }=0 \hbox{ in trace sense}\}. \label{DT} \end{equation} We define $$ D(T_i)= D_Q^r $$ as the domain of $T_i$ $(i=1,\dots ,r)$. An essential special case of (\ref{SY}) is that with polynomial nonlinearity \begin{equation} f_i(x,u_1,\dots ,u_r )= \sum_{|j|\le s_i} c^{(i)}_{j_1,\dots ,j_r}(x) u_1^{j_1}\dots u_r^{j_r} \label{Pol} \end{equation} that fulfils condition (C4), where $s_i \in\mathbb{N}^+$, $s_i \le p-1$, $c^{(i)}_{j_1,\dots ,j_r}\in C(\overline\Omega)$ and $|j|\equiv j_1+\dots j_r$ for $j=(j_1,\dots ,j_r)\in \mathbb{N}^r$. This occurs in steady states or in time discretization of autocatalytic reaction-diffusion systems. (Then $a_i$ and $c^{(i)}_{j_1,\dots ,j_r}$ are generally constant). \subsection*{(a) Energy spaces} The construction of the gradient method requires the introduction and the study of some properties of energy spaces of the Laplacian. \begin{definition} \label{def22} \rm Let $$ Bu\, \equiv \, -\Delta u + cu , $$ defined for $u\in D(B)={D_Q}$ (see (\ref{DT})), where $c={\kappa \over m}$ $(\ge 0)$ with $m$ and $\kappa$ from conditions (C3)-(C4) \, (i.e. $B=-\Delta$ except the case $\alpha\equiv 0$). \end{definition} \begin{remark} \label{rem21} \rm It can be seen in the usual way that $B$ is symmetric and strictly positive in the real Hilbert space $L^2(\Omega)$. \end{remark} \begin{corollary} \label{cor23} \begin{enumerate} \item[\rm (a)] The eigenvalues $\lambda_i$ $(i\in\mathbb{N}^+)$ of $B$ are positive. \item[\rm (b)] We have $ \int_\Omega (Bu)u \ge \lambda_1 \int_\Omega u^2 $ $(u\in {D_Q})$ where $\lambda_1>0$ is the smallest eigenvalue of $B$. \end{enumerate} \end{corollary} \begin{definition} \label{def23} \rm Denote by $H^1_Q(\Omega)$ the $energy \ space$ of $B$, i.e. $ H^1_Q(\Omega)= H_B $ (cf. Definition \ref{def21}). Due to the divergence theorem we have \begin{equation} \langle u,v \rangle _{H^1_Q} \equiv \, \langle u,v \rangle_B = \int_\Omega \bigl( \nabla u \cdot \nabla v + cuv \bigr)\, dx + \int_{\Gamma_\beta} {\alpha \over \beta}uv \, d\sigma \quad (u,v\in D). \label{UVQ} \end{equation} \end{definition} \begin{remark} \label{rem22} \rm Using Corollary \ref{cor23}, we can deduce the following properties: \begin{enumerate} \item[\rm (a)] $ (1+\lambda^{-1}) \|u\|^2_{H^1_Q} \ge \|u\|_{H^1(\Omega)} \equiv \int_\Omega \bigl( |\nabla u|^2 + u^2 \bigr)\, dx \quad (u\in H^1_Q(\Omega)).$ \item[\rm (b)] $H^1_Q(\Omega)\subset H^1(\Omega)$. \item[\rm (c)] (\ref{UVQ}) holds for all $u,v\in H^1_Q(\Omega)$. \end{enumerate} \end{remark} \begin{remark} \label{rem23} \rm Remark \ref{rem22} (b) implies that the Sobolev embedding theorem \cite{Ad} holds for $H^1_Q(\Omega)$ in the place of $H^1(\Omega)$. Namely, for any $p\ge 2$ if $N=2$, and for $2\le p\le {2N \over N-2}$ \, if $N>2$, there exists $K_{p,\Omega}>0$ such that \begin{equation} H^1_Q(\Omega) \subset L^p(\Omega), \quad \|u\|_{L^p(\Omega)} \le K_{p,\Omega} \|u\| _{H^1_Q} \quad (u\in H^1_Q(\Omega)). \label{Sob} \end{equation} \end{remark} \begin{definition} \label{def24} \rm The product spaces $L^2(\Omega)^r$ and $H^1_Q(\Omega)^r$ are endowed with the norms $$ \|u\|_{L^2(\Omega)^r}\equiv \bigl( \sum_{i=1}^r \|u_i \|^2_{L^2(\Omega)} \bigr) ^{1/2} \quad \hbox{and} \quad \|u\|_{H^1_Q(\Omega)^r}\equiv \bigl( \sum_{i=1}^r \|u_i \|^2 _{H^1_Q} \bigr) ^{1/2}, $$ respectively, where $u=(u_1,\dots ,u_r)$ and $\| \cdot \|_{H^1_Q} = \| \cdot \|_{H^1_Q(\Omega)}$ for brevity as in Def. 2.3. \end{definition} \subsection*{(b) The convergence result} \begin{theorem} \label{th23} Under the conditions (C1)-(C4) the following results hold. \\ {\rm (1)} The system (\ref{SY}) has a unique weak solution $u^\ast=(u^\ast_1, \dots ,u^\ast_r)\in H^1_Q(\Omega)^r$. \\ {\rm (2)} Let $u^0_i\in D_Q$ $(i=1,\dots ,r)$ and \begin{equation} \begin{aligned} M&= \max \{1,\eta \}m' + \kappa'\lambda_1^{-1}\\ &\quad + \gamma K_{p,\Omega}^p \mu_p \big( \|u^0\|_{H^1_Q(\Omega)^r} + m^{-1}\lambda_1^{-1/2} \|T(u^0)-g\|_{L^2(\Omega)^r} \big)^{p-2} \end{aligned}\label{Mdef} \end{equation} (with $m$, $m'$, $\eta$ from condition (C3), $p$, $\kappa'$, $\gamma$ from (C4), $K_{p,\Omega}$ from Remark \ref{rem23}, $\lambda_1$ from Corollary \ref{cor23} and \, $\mu_p=\max \{1, r^{(4-p)/2} \}$). \item[ ] Let \begin{equation} u^{k+1}_i= u^k_i- \frac{2}{M+m} z^k_i \quad (k\in\mathbb{N} , \, i=1,\dots ,r) \label{iter-a} \end{equation} where \begin{equation} g^k_i= T_i(u^k)-g_i \quad (k\in\mathbb{N} , \, i=1,\dots ,r) \label{iter-b} \end{equation} and $z^k_i\in D_Q$ is the solution of the auxiliary problem \begin{equation} \begin{gathered} (-\Delta+ c )z^k_i= g^k_i \cr \alpha(x) z^k_i +\beta(x) \partial_\nu z^k_i {}\big|_{\partial \Omega }=0\, . \end{gathered} \label{iter-c} \end{equation} (We solve Poisson equations $-\Delta z^k_i = g^k_i$\, , owing to $c=0$, except the case of the Neumann problem.) \item[ ] Then the sequence $(u^k)=(u^k_1, \dots , u^k_r)\subset D_Q^r$ converges to $u^\ast$ according to the linear estimate $$ \|u^k-u^\ast\|_{H^1_Q(\Omega)^r}\le \frac{1}{m\sqrt{\lambda_1}} \|T(u^0)-g\|_{{L^2(\Omega)}^r} \big( \frac{M-m}{M+m} \big)^k \quad (k\in\mathbb{N}^+). $$ (Owing to Remark \ref{rem22} this also means convergence in the usual $H^1(\Omega)$ norm.) \end{theorem} \begin{proof} (a) First we remark the following facts: condition (C4) implies that for all $i,k=1,\dots ,r$ and $(x,\xi )\in \overline\Omega \times \mathbb{R}^r$ $$ \left | \partial_{\xi_k} f_i (x, \xi )\right | \le \kappa'+\gamma \sum_{j=1}^r |\xi_j|^{p-2} \, .$$ Hence from Lagrange's inequality we have for all $i=1,\dots ,r$, $(x,\xi) \in \overline \Omega \times {\bf R}^r $ \begin{equation} |f_i (x,\xi )| \le | f_i(x,0)| + \Big( \kappa'+\gamma \sum_{j=1}^r |\xi_j|^{p-2}\Big) \sum_{k=1}^{r} |\xi _k| \le \kappa''+\gamma' \sum_{j=1}^r |\xi_j|^{p-1} \label{Lag} \end{equation} with suitable constants $\kappa'', \gamma '>0$. \noindent (b) To prove our theorem, we have to check conditions (i)-(iv) of Theorem 2.2 and (v)' of Corollary \ref{cor22} in our setting in the real Hilbert space $H=L^2(\Omega)$. \begin{enumerate} \item[\rm (i)] For any $u\in D_Q^r$ we have $$ |T_i(u)| \le \sum_{k=1}^{N} \left( |\partial_k a_i \, \partial_k u_i| + |a_i\, \partial^2_k u_i| \right) + |f_i(x,u_1,\dots ,u_r )|\, . $$ Here $\partial_k a_i$ and $a_i$ are in $C(\overline\Omega)$, $\partial_k u_i$ and $\partial^2_k u_i$ are in $L^2(\Omega)$, hence the sum term is in $L^2(\Omega)$. Further, assumption (C4) implies $2p-2< {2N \over N-4}$, hence $H^2(\Omega) \subset L^{2p-2}(\Omega)$ \, \cite{Ad}. Thus (\ref{Lag}) yields $|f_i(x,u_1,\dots ,u_r )|\le \kappa''+\gamma' \sum_{j=1}^r |u_j|^{p-1} \, \in L^2(\Omega).$ That is, $T_i$ maps indeed into $L^2(\Omega)$. Further, assumption s (C1)-(C2) imply that for any $g\in L^2(\Omega)$ the weak solution of $-\Delta z +cz =g$ \, with $\alpha z +\beta \partial_\nu z \big|_{\partial \Omega }=0$ \, is in $H^2(\Omega)$ \, \cite{Enc}, i.e. $R(B)=L^2(\Omega)$. Hence $R(B)\supset R(T_i)$ holds. \item[\rm (ii)] For any $u\in D_Q^r$, $v\in D_Q$ and $i=1,\dots ,r$ the divergence theorem yields (similarly to (\ref{UVQ})) \begin{equation} \begin{aligned} \langle B^{-1}T_i(u),v\rangle_{H^1_Q} &= \int_\Omega T_i(u)v\\ &= \int_\Omega \bigl( a_i \, \nabla u_i \cdot \nabla v + f_i(x,u)v \bigr)\, dx + \int_{\Gamma_\beta} a_i {\alpha \over \beta}u_iv \, d\sigma \, . \end{aligned}\label{obtint} \end{equation} Let us put arbitrary $u\in H^1_Q(\Omega)^r, v\in H^1_Q(\Omega)$ in (\ref{obtint}). Setting $Q(u)\equiv \gamma' \sum_{j=1}^r |u_j|^{p-1}= \gamma' \sum_{j=1}^r |u_j|^{p/q}$, where $p^{-1}+q^{-1}=1$, we have $|f_i(x,u)|\le \kappa''+Q(u)$ from (\ref{Lag}) where $Q(u)\in L^q(\Omega)$. Then (\ref{obtint}) can be estimated by the expression \begin{gather*} \max_{\overline\Omega} a_i \, \|\nabla u \|_{L^2(\Omega)} \|\nabla v \|_{L^2(\Omega)} + \kappa''|\Omega|^{1/2} \|v \|_{L^2(\Omega)} \\ + \|Q(u) \|_{L^q(\Omega)} \|v \|_{L^p(\Omega)} + \eta \max_{\Gamma_\beta} a_i\, \|u\|_{L^2(\partial\Omega)} \|v\|_{L^2(\partial\Omega)}, \end{gather*} where $|\Omega|$ is the measure of $\Omega$. Using (\ref{Sob}) and the continuity of the trace operator, we see that for any fixed $u\in H^1_Q(\Omega)^r$ the expression (\ref{obtint}) defines a bounded linear functional on $H^1_Q(\Omega)^r$. Hence (using Riesz's theorem) we define the operator $F_i: H^1_Q(\Omega)^r \to H^1_Q(\Omega)$ by the formula $$ \langle F_i(u),v\rangle _{H^1_Q} = \int_\Omega \bigl( a_i \, \nabla u_i \cdot \nabla v + f_i(x,u)v \bigr)\, dx + \int_{\Gamma_\beta} a_i {\alpha \over \beta}u_iv \, d\sigma\,, $$ $u\in H^1_Q(\Omega)^r, v\in H^1_Q(\Omega)$. \end{enumerate} For $u\in H^1_Q(\Omega)^r$ let $S_i(u) \in B\left( H^1_Q(\Omega)^r, H^1_Q(\Omega)\right )$ be the bounded linear operator defined by \begin{equation} \langle S_i(u)h,v\rangle_{H^1_Q} = \int_\Omega \bigl( a_i \, \nabla h_i \cdot \nabla v + \sum_{k=1}^r \partial_{\xi_k} f_i (x, u) h_k v \bigr)\, dx + \int_{\Gamma_\beta} a_i {\alpha \over \beta}h_iv \, d\sigma \label{SIUHV} \end{equation} $u\in H^1_Q(\Omega)^r, v\in H^1_Q(\Omega)$. The existence of $S_i(u)$ is provided by Riesz's theorem similarly as above. Now having the estimate \begin{align*} & \int_\Omega |\partial_{\xi_k} f_i (x, u) h_k v|\, dx \\ &\le \kappa' \|h_k \|_{L^2(\Omega)} \|v \|_{L^2(\Omega)} + \gamma \big\| \sum_{k=1}^r |u_j|^{p-2} \big\| _{L^{p \over p-2}(\Omega)} \|h_k \|_{L^p(\Omega)} \|v \|_{L^p(\Omega)} \end{align*} for the terms with $f_i$, using $( {p \over p-2})^{-1}+p^{-1}+p^{-1}=1$. We will prove that $F_i$ is G\^ateaux differentiable $(i=1,\dots ,r)$, namely, $$ F_i ' (u)=S_i(u) \quad \left ( u\in H^1_Q(\Omega)^r \right )\, . $$ Let $u,h \in H^1_Q(\Omega)^r$, further, $\mathcal{E}\equiv \{ v\in H^1_Q(\Omega): \|v\|_{H^1_Q(\Omega)} = 1\}$ and \begin{align*} \delta^i_{u,h}(t)&\equiv \frac{1}{t}\| F_i(u+th)-F_i(u)- tS_i(u)h\|_{H^1_Q(\Omega)} \\ &= \frac{1}{t} \sup_{v\in \mathcal{E}} \langle F_i(u+th)- F_i(u)-tS_i(u)h,v\rangle _{H^1_Q(\Omega)}\, . \end{align*} Then, using linearity, we have \begin{align*} \delta^i_{u,h}(t) & =\frac{1}{t}\sup_{v\in \mathcal{E}} \int_\Omega \Big( f_i(x,u+th)-f_i(x,u)- t \sum_{k=1}^r \partial_{\xi_k} f_i (x, u) h_k\Big) v\, dx\\ &= \sup_{v\in \mathcal{E}} \int_\Omega \sum_{k=1}^r \left( \partial_{\xi_k} f_i (x,u+\theta th) - \partial_{\xi_k} f_i (x, u) \right) h_k v\, dx\\ &\le \sup_{v\in \mathcal{E}} \sum_{k=1}^r \Big( \int_\Omega |\partial_{\xi_k} f_i (x,u+\theta th) - \partial_{\xi_k} f_i (x, u)|^{p \over p-2} dx \Big) ^{p-2 \over p}\\ &\quad\times \|h_k \|_{L^p(\Omega)} \|v \|_{L^p(\Omega)} \, . \end{align*} Here $\|v \|_{L^p(\Omega)} \le K_{p,\Omega}$ from (\ref{Sob}), further, $|\theta th| \le |th| \to 0$ a.e. on $\Omega $, hence the continuity of $\partial _{\xi _k}f_i$ implies that the integrands converge a.e. to 0 when $t\to 0$. The integrands are majorized for $t\le t_0$ by $$ \Big| 2(\kappa' + \gamma \sum_{j=1}^r |u_j+t_0 h_j|^{p-2}) \Big| ^{p \over p-2} \le \tilde \kappa + \tilde \gamma \sum_{j=1}^r |u_j+t_0 h_j|^p $$ in $L^1(\Omega )$, hence, by Lebesgue's theorem, the obtained expression tends to 0 when $t \to 0$. Thus $$ \lim_{t\rightarrow 0} \delta^i_{u,h}(t) =0\, . $$ \begin{enumerate} \item[\rm (iii)] It is proved similarly to (ii) that for fixed functions $ u, k,w\in H^1_Q(\Omega)^r, h \in H^1_Q(\Omega)$ the mapping $s,t \mapsto F_i'(u+sk+tw)h$ is continuous from ${\bf R}^2$ to $H^1_Q(\Omega)$. Namely, \begin{align*} \omega _{u,k,w,h}(s,t)&\equiv \| F_i'(u+sk+tw)h- F_i'(u)h\|_{H^1_Q(\Omega)} \\ &= \sup_{v\in \mathcal{E}} \langle F_i'(u+sk+tw)h- F_i'(u)h,v\rangle _{H^1_Q(\Omega)}\\ &= \sup_{v\in \mathcal{E}} \int_\Omega \sum_{k=1}^r \left( \partial_{\xi_k} f_i (x,u+sk+tw) - \partial_{\xi_k} f_i (x, u) \right) h_k v\, dx \, . \end{align*} Then we obtain just as above (from the continuity of $\partial_{\xi_k} f_i$ and Lebesgue's theorem) that $$\lim_{s,t\rightarrow 0} \omega _{u,k,w,h} (s,t) =0\, .$$ \item[\rm (iv)] It follows from $F_i'(u)=S_i(u)$ in (\ref{SIUHV}) and from the assumed symmetry of the Jacobians $\partial_\xi f (x,\xi)$ that for any $u,h,v\in H^1_Q(\Omega)^r$ $$ \sum_{i=1}^{r} \langle F_i'(u)h,v_i\rangle _{H^1_Q(\Omega)} = \sum_{i=1}^{r} \langle h_i,F_i'(u)v \rangle _{H^1_Q(\Omega)}. $$ \item[\rm (v)] For any $u,h \in H^1_Q(\Omega)^r$ we have \begin{align*} & \sum_{i=1}^{r} \langle F'_i(u)h,h_i\rangle _{H^1_Q(\Omega)} \\ &= \int_\Omega \bigl( \sum_{i=1}^{r} a_i |\nabla h_i|^2 + \sum_{i,k=1}^r \partial_{\xi_k} f_i (x, u) h_k h_i \bigr)\, dx + \int_{\Gamma_\beta} {\alpha \over \beta} \sum_{i=1}^{r} a_i h_i^2 \, d\sigma \, . \end{align*} \end{enumerate} Hence from assumptions (C3)-(C4) we have \begin{align*} &\sum_{i=1}^{r} \langle F'_i(u)h,h_i\rangle _{H^1_Q(\Omega)} \\ &\ge m \int_\Omega \sum_{i=1}^{r} |\nabla h_i|^2 \, dx + \kappa \int_\Omega \sum_{i=1}^{r} h_i^2 \, dx + m \int_{\Gamma_\beta} {\alpha \over \beta} \sum_{i=1}^{r} h_i^2 \, d\sigma \\ &= m\|h\|_{H^1_Q(\Omega)^r} ^2 \end{align*} using $\kappa=cm$ (see Def.2.2). Further, \begin{equation} \begin{aligned} \sum_{i=1}^{r} \langle F'_i(u)h,h_i\rangle _{H^1_Q(\Omega)} &\le m' \sum_{i=1}^{r} \int_\Omega |\nabla h_i|^2 \, dx + \eta m'\sum_{i=1}^{r} \int_{\Gamma_\beta} h_i^2 \, d\sigma \\ &\quad + \int_\Omega \Big( \kappa'+\gamma \sum_{j=1}^r |u_j|^{p-2}\Big) \sum_{k=1}^r h_k^2\, dx \, . \end{aligned}\label{uhp1} \end{equation} Here $$ \kappa' \sum_{k=1}^r \int_\Omega h_k^2\, dx \le {\kappa' \over \lambda_1} \|h\|^2_{H^1_Q(\Omega)^r} $$ from Corollary \ref{cor23} (b). Further, from H\"older's inequality, using ${p-2\over p}+ {2\over p}=1$, we obtain \begin{align*} \sum_{j,k=1}^r \int_\Omega |u_j|^{p-2} h_k^2\, dx &\le \sum_{j,k=1}^r \Big[ \int_\Omega \left(|u_j|^{p-2}\right)^{p\over p-2}\Big]^{p-2\over p} \Big[\int_\Omega \left(h_k^2\right)^{p/ 2}\Big]^{2/ p}\\ &=\sum_{j,k=1}^r \|u_j\|^{p-2}_{L^p(\Omega)} \|h_k\|_{L^p(\Omega)}^2\\ &= \Big( \sum_{j=1}^r \|u_j\|^{p-2}_{L^p(\Omega)} \Big) \Big( \sum_{k=1}^r \|h_k\|_{L^p(\Omega)}^2 \Big)\, . \end{align*} An elementary extreme value calculation shows that for $x\in\mathbb{R}^r$, $\sum_{j=1}^r x_j^2 = R^2 $ the values of $ \left( \sum_{j=1}^r |x_j|^{p-2} \right) ^{2 \over p-2} $ lie between $R^2$ and $r^{4-p \over p-2}R^2$, i.e. $$ \sum_{j=1}^r |x_j|^{p-2} \le \mu_p \Big( \sum_{j=1}^r |x_j|^2 \Big) ^{p-2 \over 2} $$ where $\mu_p = \max \{ 1, r^{4-p \over p-2} \} $. Hence \begin{align*} \sum_{j,k=1}^r \int_\Omega |u_j|^{p-2} h_k^2\, dx &\le \mu_p \Big( \sum_{j=1}^r \|u_j\|^2_{L^p(\Omega)} \Big)^{p-2 \over 2} \Big( \sum_{k=1}^r \|h_k\|_{L^p(\Omega)}^2 \Big)\\ &\le \mu_p K_{p,\Omega}^p \Big( \sum_{j=1}^r \|u_j\|^2_{H^1_Q} \Big)^{p-2 \over 2} \Big( \sum_{k=1}^r \|h_k\|_{H^1_Q}^2 \Big)\\ &= \mu_p K_{p,\Omega}^p \|u\|_{H^1_Q(\Omega)^r}^{p-2} \|h\|_{H^1_Q(\Omega)^r}^2\, . \end{align*} Summing up, (\ref{uhp1}) yields $$ \sum_{i=1}^{r} \langle F'_i(u)h,h_i\rangle _{H^1_Q(\Omega)} \le\, M(u) \|h\|^2_{H^1_Q(\Omega)^r} \quad (u,h\in H^1_Q(\Omega)^r) $$ with $$ M(u) = \max \{1,\eta \}m' + \kappa'\lambda_1^{-1} + \gamma K^p_p(\Omega) \mu_p \|u\|_{H^1_Q(\Omega)^r}^{p-2}\, . $$ Since Corollary \ref{cor23} (b) implies $\|u\|_{H^1_Q} \le \lambda_1^{-1/2} \|Bu\|_{L^2(\Omega)}$ $(u\in D_Q)$, therefore the radius $r_0= m^{-1} \|F(u^0)-b\|_{H^1_Q(\Omega)^r} $ (defined in Corollary \ref{cor22}) fulfils $$ r_0 \le m^{-1}\lambda_1^{-1/2} \bigl( \sum_{i=1}^r \|T_i(u^0)-g_i\|_{L^2(\Omega)}^2 \bigr)^{1/2}\, , $$ using $BF_{i \mid D} = T_i$. Hence for $u\in B(u^0,r_0)=\{u\in H^1_Q(\Omega)^r:\, \|u-u^0\|_{H^1_Q(\Omega)^r} \le r_0 \}$ we have $$ \sum_{i=1}^{r} \langle F'_i(u)h,h_i\rangle _{H^1_Q(\Omega)} \le\, M \|h\|^2_{H^1_Q(\Omega)^r} \quad (u\in B(u^0,r_0), h\in H^1_Q(\Omega)^r) $$ with $M$ defined in (\ref{Mdef}). \end{proof} \begin{remark} \label{rem24} \rm Theorem \ref{th23} holds similarly in the following cases: \\ (a)\; with other smoothness assumption s on $\partial\Omega$ and the coefficients of $T$, when the inclusion $R(T_i)\subset R(B)$ is fulfilled with suitable domain $D(T_i)$ of $T_i$ instead of $D_Q^r$ (cf. (\ref{DT})). \\ (b)\; with more general linear part $ -\mathop{\rm div} (A_i(x)\nabla u_i)$ of $T_i$, where $A_i\in C^1(\overline\Omega, \mathbb{R}^{N\times N})$, in the case of Dirichlet boundary condition. \end{remark} The above theorem is the extension of the cited earlier results on the gradient method to system (\ref{SY}). We note that the {\it conjugate gradient method } might be similarly extended to (\ref{SY}), following its application in \cite{cgmn} for a single Dirichlet problem. As mentioned earlier, the CGM constructs two simultaneous sequences, and it improves the convergence ratio of the gradient method to $(\sqrt{M}-\sqrt{m})(\sqrt{M}+\sqrt{m})$ at the price of an %on the other hand, the extra work which comprises the calculation of required integrals and numerical root finding for the stepsize. Compared to {\it Newton-like methods} (which can provide higher order convergence than linear), we emphasize that in the iteration of Theorem 2.3 the auxiliary problems are of fixed (Poisson) type, whereas Newton-like methods involve stepwise different linearized problems with variable coefficients. Hence in our iteration one can exploit the efficient solution methods that exist for the Poisson equation. (More discussion on this will be given in subsection \ref{seceffpo}.) \subsection{Sobolev embedding properties} The construction of the sequence (\ref{iter-a}) in Theorem 2.3 needs an explicit value of the constant $M$ in (\ref{Mdef}). The parameters involved in $M$ are defined in conditions (C3)-(C4) only with the exception of the embedding constants $K_{p,\Omega}$ in (\ref{Sob}). (The eigenvalue $\lambda_1$ fulfils $\lambda_1= K_{2,\Omega}^{-2}$ by virtue of Corollary \ref{cor23} (b).) Consequently, in order to define the value of $M$, the exact value or at least a suitable estimate is required for the constants $K_{p,\Omega}$. Although most of the exact constant problems have been solved in $\mathbb{R}^n$, even in the critical exponent case (see \cite{Ros, Tal}), for bounded domains there is not yet complete answer. For $n\ge 3$ and for small square in the case $n=2$, the embedding constants of $H^1(\Omega)$ to $L^p(\Omega)$ are given in \cite{Bur1, Bur2}; an estimate is given for functions partly vanishing on the boundary for the critical exponent case in \cite{Li}. Consequently, a brief study of the embeddings is worth wile to obtain estimates of the embedding constants which are valid for $n=2$. Our estimates, presented in two dimensions, take into account the boundary values of the functions. Besides the constants $K_{p,\Omega}$ in (\ref{Sob}), for any set $\Gamma\subset \partial\Omega$ we denote by $K_{p,\Gamma}$ the embedding constant in the estimate $$ \| u\big|_{\Gamma} \|_{L^p(\Gamma)} \le K_{p,\Gamma} \|u\|_{H^1_Q} \quad (u\in H^1_Q(\Omega)). $$ \begin{lemma} \label{lem21} Let $I=[a,b]\times[c,d]\subset\mathbb{R}^2$, $p_i\ge 1$ \ $(i=1,2)$. The boundary $\partial I$ is decomposed into $\Gamma_1 = \{ a,b \}\times [c,d]$ and $\Gamma_2 = [ a,b ]\times \{c,d\}$. Then $$ K_{p_1+p_2,I}^{p_1+p_2} \le \frac{1}{2} \big( K_{p_1,\Gamma_1}^{p_1} + p_1 K_{2(p_1-1),I}^{p_1-1} \big) \big( K_{p_2,\Gamma_2}^{p_2} + p_2 K_{2(p_2-1),I}^{p_2-1} \big)\, . $$ \end{lemma} \begin{proof} Let $u\in H^1(I)$. We define the function $u_a(y)\equiv u(a,y)$, and similarly $u_b$, $u_c$ and $u_d$. Then for any $x,y\in I$ we have \begin{equation} \begin{aligned} |u(x,y)|^{p_1} &= |u_a(y)|^{p_1} + p_1 \int_a^x |u(\xi,y)|^{p_1-2} u(\xi,y)\partial_1 u(\xi,y)\, d\xi \\ &\le |u_a(y)|^{p_1} + p_1 \int_a^b |u|^{p_1-1} |\partial_1 u| \, dx\, . \end{aligned}\label{up1} \end{equation} Similarly, we obtain \begin{equation} |u(x,y)|^{p_2} \le |u_c(y)|^{p_2} + p_2 \int_c^d |u|^{p_2-1} |\partial_2 u| \, dy\, . \label{up2} \end{equation} Multiplying (\ref{up1}) and (\ref{up2}) and then integrating over $I$, we obtain \begin{align*} \int_I |u|^{p_1+p_2} &\le \Big( \int_c^d |u_a|^{p_1} + p_1 \int_I |u|^{p_1-1} |\partial_1 u| \Big) \Big( \int_a^b |u_c|^{p_2} +p_2 \int_I |u|^{p_2-1} |\partial_2 u| \Big)\\ &\le \Big( \int_c^d |u_a|^{p_1} + p_1 \|u\|^{p_1-1}_{L^{2(p_1-1)}(I)} \|\partial_1 u\|_{L^2(I)} \Big)\\ &\quad\times \Big( \int_a^b |u_c|^{p_2} + p_2 \|u\|^{p_2-1}_{L^{2(p_2-1)}(I)} \|\partial_2 u\|_{L^2(I)} \Big)\, . \end{align*} The same holds with $u_b$ and $u_d$ instead of $u_a$ and $u_b$. Using the elementary inequality \begin{align*} &(\alpha_1 + r_1 \gamma_1)(\alpha_2 + r_2 \gamma_2) + (\beta_1 + r_1 \gamma_1)(\beta_2 + r_2 \gamma_2)\\ &\le (\alpha_1 + \beta_1 + r_1 \sqrt{\gamma_1^2 + \gamma_2^2} ) (\alpha_2 + \beta_2 + r_2 \sqrt{\gamma_1^2 + \gamma_2^2} ) \end{align*} for $\alpha_i,\beta_i,r_i,\gamma_i\ge 0$ $(i=1,2)$, we obtain \begin{align*} 2 \int_I |u|^{p_1+p_2} & \le \Big( \int_c^d (|u_a|^{p_1} + |u_b|^{p_1}) + p_1 \|u\|^{p_1-1}_{L^{2(p_1-1)}(I)} \|\nabla u\|_{L^2(I)} \Big)\\ &\quad\times \Big( \int_a^b (|u_c|^{p_2} +|u_d|^{p_2}) + p_2 \|u\|^{p_2-1}_{L^{2(p_2-1)}(I)} \|\nabla u\|_{L^2(I)} \Big)\\ & \le \left( K_{p_1,\Gamma_1}^{p_1} + p_1 K_{2(p_1-1),I}^{p_1-1} \right) \left( K_{p_2,\Gamma_2}^{p_2} + p_2 K_{2(p_2-1),I}^{p_2-1} \right) \|u\|_{H^1_Q}^{p_1+p_2}\, . \end{align*} \end{proof} \begin{corollary} \label{cor24} Let $\Omega\subset\mathbb{R}^2$ with $\partial\Omega\in C^1$ and let us consider Dirichlet boundary condition s in {\rm (\ref{SY})}, i.e. $Qu\equiv u$ and $H^1_Q(\Omega) =H^1_0(\Omega)$. Then \begin{equation} K_{p_1+p_2,\Omega}^{p_1+p_2} \le \frac{p_1 p_2}{2} K_{2(p_1-1),\Omega}^{p_1-1} K_{2(p_2-1),\Omega}^{p_2-1} \, . \label{Cor24} \end{equation} \end{corollary} \begin{proof} $\Omega$ is included in some $I=[a,b]\times[c,d]$, and for any $u\in H^1_0(\Omega)$ its extension $\tilde u\in H^1_0(I)$ is defined as zero on $I\setminus \Omega$. Then for any $p\ge 1$ we have $K_{p,\Gamma_i}=0$ and $K_{p,\Omega} = K_{p,I}$. \end{proof} The case when $p$ is an even integer is of particular importance since we have this situation in the case of (\ref{Pol}) owing to (C4). Then the functional inequality (\ref{Cor24}) leads directly to an estimate: \begin{corollary} \label{cor25} Let $\lambda_1>0$ be the smallest eigenvalue of $-\Delta$ on $H^1_0(\Omega)$. Then \begin{enumerate} \item[\rm (a)] $K_{2n,\Omega} \le 2^{-1/2} \left( {2 \over \lambda_1} \right)^{1/2n} (n!)^{1/n} \quad (n\in\mathbb{N}^+)$; \item[\rm (b)] $K_{2n,\Omega} \le 0.63 b_n n \quad (n\in\mathbb{N}^+, n\ge 2)$ where $b_n=(2/ \lambda_1)^{1/2n}$ (and thus $\lim b_n=1$). \end{enumerate} \end{corollary} \begin{proof} (a) Let $h(p)=K_{p,\Omega}^p$ $(p\ge 1)$. Then for $n\in\mathbb{N}^+$ (\ref{Cor24}) implies the recurrence $$ h(2n)\le h(2n-2){n^2 \over 2}, $$ hence $$ h(2n) \le h(2) \frac{ (n!)^2 }{ 2^{n-1} } = {2 \over \lambda_1} \frac{ (n!)^2 }{2^n}\, , $$ since Corollary \ref{cor23} gives $K_{2,\Omega}= \lambda_1^{-1/2}$. \noindent (b) The estimate $(n!)^{1/n}\le 0.891\, n$ $(n\ge 2)$ is used. \end{proof} The boundary embedding constants $K_{p,\Gamma_i}$ can be estimated in terms of suitable $K_{p'}(I)$ as follows. \begin{lemma} \label{lem22} Let $I$ and $\Gamma_i$ $(i=1,2)$ be as in Lemma \ref{lem21}, $p\ge 1$. Then $$ K_{p,\Gamma_i}^p \le {2 \over b-a} K_{p,I}^p + p\sqrt{2} K_{2(p-1),I}^{p-1} \, . $$ \end{lemma} \begin{proof} We prove the lemma for $\Gamma_1$. Similarly to Lemma \ref{lem21} we have $$ |u_a(y)|^p \le \,|u(x,y)|^p + p \int_a^b |u|^{p-1} |\partial_1 u| \, dx \quad (x,y\in I) . $$ Integrating over $I$, we obtain \begin{align*} (b-a) \int_c^d |u_a|^p &\le\int_I |u|^p + p(b-a) \int_I |u|^{p-1} |\partial_1 u| \\ &\le\|u\|_{L^p(\Omega)}^p + p(b-a) \|u\|_{L^{2(p-1)}(\Omega)}^{p-1} \|\partial_1 u\|_{L^2(\Omega)}\, , \end{align*} and similarly for $u_b$. Hence, summing up and using \, $ \|\partial_1 u\|_{L^2(\Omega)} + \|\partial_2 u\|_{L^2(\Omega)} \le \sqrt{2} \|\nabla u\|_{L^2(\Omega)} \le \sqrt{2} \|u\|_{H^1_Q} $, we have $$ \|u\|_{L^p(\Gamma_1)}^p \le \big( {2 \over b-a} K_{p,I}^p + p\sqrt{2} K_{2(p-1),I}^{p-1} \big) \|u\|_{H^1_Q}^p\, . $$ \end{proof} \begin{corollary} \label{cor26} For any $n\in\mathbb{N}^+$ we have $$ K_{2n,\Gamma_i} \le 1.26 c_n n\, , $$ where $c_n = {1 \over 2} \left( {4 \over {\lambda_1 (b-a)} } + {4^{n+1} \over \sqrt{2\lambda_1} } \right)^{1/2n}$ (and thus $\lim c_n=1$). \end{corollary} The proof of this corollary follow using Corollary \ref{cor25} (a) and again $n!< (0.891 n)^n$. \begin{remark} \label{rem25} \rm Lemmas 2.1 and 2.2 may be extended from the interval case to other domains, depending on the actual shape of $\Omega$, if portions $\Gamma$ of $\partial\Omega$ are parametrized as e.g. $t\mapsto (t,\varphi(t))$ and inequalities of the type $ \int_\alpha^\beta |u|^p\, dx\le \int_\alpha^\beta |u(t,\varphi(t))|^p (1+\varphi'(t))^{1/2}\, dt= \int_\Gamma |u|^p$ are used. Then estimates can be derived depending on the portions $\Gamma_i$ of $\partial\Omega$ where $u\big|_{\Gamma_i} =0$ (i.e. $K_{p,\Gamma_i}=0$). The detailed investigation is out of the scope of this paper. (For a model problem a calculation of the corresponding estimate for mixed boundary condition s will be given in section \ref{secex}.) \end{remark} \section{Implementation of the method} In Section \ref{secgen}, the Sobolev space gradient method was developed for systems of the form (\ref{SY}). Thereby a theoretical iteration is executed directly for the original boundary-value problem in the Sobolev space, and it is shown that the iteration converges in the corresponding energy norm. One of the main features of this approach is the reduction of computational questions to those arising for the auxiliary linear Poisson problems. Namely, in the application to a given system of boundary-value problems one's task is to choose a numerical method %suitable(FDM, FEM) for the Poisson problems and solve the latter to some suitable accuracy. This means that from now on two issues have to be handled: from the aspect of convergence, error control for the stepwise solution of the auxiliary problems, and from the aspect of cost, the efficient solution of the Poisson equations. Section \ref{secpreli} is devoted to these topics. First a discussion of corresponding error estimates is given for the numerically constructed sequences. Then we will refer very briefly to some efficient Poisson solvers. Here we note that the efficiency of the whole iteration much relies on the fact that all the linear problems are of the same Poisson type, for which efficient solvers are available. In Sections \ref{secdirect}--\ref{secex} we consider the simplest case of realization as an elementary illustration of the theoretical results. This suits the semilinear setting of this paper and involves the case of polynomial nonlinearity (\ref{Pol}), connected to reaction-diffusion systems. Namely, on some special domains the GM is applied in effect directly to the BVP itself since the Poisson equations are solved exactly. This is due to keeping the iteration in special function classes. We give a brief summary of some cases of such special domains. Then the paper is closed with an example that illustrates the convergence result. \section{Error control and efficiency} \label{secpreli} \subsection{Error estimates for the numerical iterations} % 3.1 \label{secerrest} The theoretical iteration $(u^k)=(u^k_1, \dots , u^k_r)$ $(k\in\mathbb{N})$, defined in Theorem \ref{th21}, can be written as \begin{gather*} u^{k+1}= u^k- \frac{2}{M+m} z^k \\ \hbox{where } z^k= \mathcal{B}^{-1}(T(u^k)-g), \end{gather*} using the notation $$ \mathcal{B}:D_Q^r\to L^2(\Omega)^r , \quad \mathcal{B}(w_1, \dots ,w_r)\equiv (Bw_1, \dots ,Bw_r) . $$ Recall that $B$ is defined on $D_Q$ (see (\ref{DT})), containing the boundary conditions, and we have $Bu\equiv -\Delta u$ if $\alpha \not\equiv 0$ and $Bu\equiv -\Delta u + cu$ if $\alpha \equiv 0$ (i.e. for Neumann BC). Any kind of numerical implementation of the GM in the Sobolev space $H^1_Q(\Omega)^r$ defines a sequence $(\overline u^k)$, constructed as follows: \begin{equation} \label{numiter} \begin{gathered} \overline{u}^0=u^0\in D; \\ \hbox{for } k\in\mathbb{N}: \quad \overline{u}^{k+1} = \overline{u}^k - \frac{2}{M+m}\overline{z}^k\, , \\ \hbox{where } \overline{z}^k \approx z^k_\ast \equiv \mathcal{B}^{-1}(T(\overline{u}^k)-g) \\ \hbox{such that } \|\overline z^k-z^k_\ast \|_{H^1_Q(\Omega)^r} \le \delta_k \end{gathered} \end{equation} where $(\delta_k)\subset\mathbb{R}^+$ is a real sequence. Then our task is to estimate $\|\overline u^k-u^\ast\|_{H^1_Q(\Omega)^r}$ in terms of the sequence $(\delta_k)$, where $u^\ast=(u^\ast_1, \dots ,u^\ast_r)\in H^1_Q(\Omega)^r$ is the weak solution of the system (\ref{SY}). We define $$ E_k\equiv\, \|u^k-\overline{u}^k \|_{H^1_Q(\Omega)^r}\, . $$ By Theorem \ref{th21} we have $$ \|\overline u^k-u^\ast\|_{H^1_Q(\Omega)^r} \le E_k + \frac{R_0}{m\sqrt{\lambda_1}} \left ( \frac{M-m}{M+m} \right )^k \quad (k\in\mathbb{N}^+) $$ where $R_0=\|T(u^0)-g\|_{L^2(\Omega)^r}$ denotes the initial residual. Hence the required error estimates depend on the behaviour of $(E_k)$. We have proved the following two results in \cite{GFEM} for a single Dirichlet problem. Since they are entirely based on the bounds $m$ and $M$ of the generalized differential operator (which is also the background for our Theorem \ref{th21}), they can be immediately formulated in our setting, too. \begin{proposition}[\cite{GFEM}] \label{prop41} For all $k\in\mathbb{N}$ $$ E_{k+1}\le \frac{M-m}{M+m} E_k + \frac{2}{M+m} \delta_k . $$ \end{proposition} \begin{corollary}[\cite{GFEM}] \label{cor41} Let $00$ be fixed, $\delta_k\le c_1 q^k$ $(k\in\mathbb{N})$. Then the following estimates hold for all $k\in\mathbb{N}^+$: \begin{enumerate} \item[(a)] if $ q> \frac{M-m}{M+m}$ then $\|\overline u^k-u^\ast\|_{H^1_Q(\Omega)^r} \le c_2 q^k$; \item[(b)] if $ q< \frac{M-m}{M+m}$ then $ \|\overline u^k-u^\ast\|_{H^1_Q(\Omega)^r} \le c_3 \left(\frac{M-m}{M+m}\right)^k$ \end{enumerate} where $c_2= {\alpha c_1 \over q-Q} + \frac{R_0}{m\sqrt{\lambda_1}}$, $c_3= {\alpha c_1 \over Q-q} + \frac{R_0}{m\sqrt{\lambda_1}}$, $Q=\frac{M-m}{M+m}$. \end{corollary} Besides these convergence results, it is also of practical interest to ensure that we arrive only in a prescribed neighbourhood of the solution. \begin{proposition} \label{prop42} Let $\varepsilon >0$ be fixed, let $\delta_k\equiv m\varepsilon $. Then we have for all $k\in\mathbb{N}^+$ \begin{equation} E_k\, < \, \varepsilon\, . \label{eke} \end{equation} \end{proposition} \begin{proof} By definition $E_0= \|u^0-\overline u^0\|_{H^1_Q(\Omega)^r}=0$. Assume that (\ref{eke}) holds for fixed $k\in\mathbb{N}^+$. Then Proposition \ref{prop41} yields $$ E_{k+1}\, < \, \frac{M-m}{M+m} \varepsilon + \frac{2}{M+m} m\varepsilon= \varepsilon \, . $$ \end{proof} \begin{corollary} \label{cor42} If $(\delta_k)$ is chosen as in Proposition \ref{prop42}, then for $(k\in\mathbb{N}^+)$, $$ \|\overline u^k-u^\ast\|_{H^1_Q(\Omega)^r} \le \varepsilon + \frac{R_0}{m\sqrt{\lambda_1}} \Big( \frac{M-m}{M+m} \Big)^k. $$ \end{corollary} \subsection{Discretization and efficient Poisson solvers} \label{seceffpo} The numerical solution of the Poisson equations is generally achieved by using some discretization, which is a finite difference method or a finite element method. In this respect we note two facts concerning the whole iteration. First, if we use one and the same fixed grid for each linear problem, then (\ref{numiter}) is equivalent to a suitably precondition ed nonlinear FEM iteration such that the precondition er is the discrete Laplacian corresponding to the grid. On the other hand, the use of variable grids provide a multilevel type iteration. We emphasize that for such iterations the convergence ratio in Theorem \ref{th21} yields an asymptotic value for those of the discretized problems, hence the numerical iterations using different grids (or grid sequences) exhibit mesh uniform convergence. The numerical efficiency of the iteration (\ref{numiter}) relies on the fact that all the linear problems are of the same Poisson type. Namely, various fast Poisson solvers are available that have been developed in the past decades. Many of these solvers were originally developed on rectangular domains and later extended to other domains via the fictitious domain approach. The fast direct solvers include the method of cyclic reduction, the fast Fourier transform and the FACR algorithm: comprehensive summaries on the solution of the Poisson equation with these methods are found in \cite{Dorr, Swarz77}. Another family of fast solvers which has undergone recent development is formed by spectral methods \cite{Funaro92}. The parallel implementation of these algorithms is also feasible. For the fictitious domain approach for general domains see \cite{BW}. \section{Direct realization on special domains} % 3.2 \label{secdirect} We now consider the system (\ref{SY}) with polynomial nonlinearity (\ref{Pol}). That is, we investigate the system \begin{equation} \label{ujSY} \begin{gathered} T_i(u_1,\dots ,u_r)\equiv\, -\mathop{\rm div} (a_i(x)\nabla u_i) + \sum_{|j|\le s_i} c^{(i)}_j(x) u_1^{j_1}\dots u_r^{j_r} = g_i(x) \\ Qu_i \equiv (\alpha(x) u_i +\beta(x) \partial_\nu u_i) {}\big|_{\partial \Omega }=0 \end{gathered} \end{equation} that fulfils conditions (C1)-(C4) (given in Section 2), where the moving index is $j=(j_1,\dots ,j_r)\in\mathbb{N}^r$ with $|j|=j_1+\dots j_r$ and we have $c^{(i)}_j \in C(\overline\Omega)$, $s_i\in \mathbb{N}^+$, $s_i\le p-1$. (This type of equations occurs e.g. in steady states or in time discretization of reaction-diffusion systems. Condition (C4) expresses that the reaction is of autocatalytic type.) The main idea is the following: on special domains, first approximating the coefficients and right-sides of (\ref{ujSY}) by appropriate algebraic /trigonometric polynomials, the iteration (\ref{iter-a}) can also be kept in a suitable class of algebraic /trigonometric polynomials. Hence the solution of the auxiliary Poisson equations can be achieved directly by executing a linear combination (or, in the case of a ball, by solving a simply structured linear system of algebraic equations). The solution of the approximate system (in which the coefficients and right-sides are approximated by polynomials) remains appropriately close to the solution of (\ref{ujSY}), depending on the accuracy of approximation. This can be formulated immediately in the case when the coefficients $a_i$ and $c_j^i$ are unchanged. (We typically have this situation with the operators $T_i$ occurring in reaction-diffusion systems, wherein $a_i$ and $c_j^i$ are generally constant.) The evident but lengthy calculation for the variable coefficient case is left to the reader. \begin{proposition} \label{prop51} Let $\varepsilon>0$. Let $\delta_i>0$ $(i=1,\dots ,r)$ such that $(\sum_{i=1}^r \delta_i^2)^{1/2}< \lambda_1^{1/2}m\varepsilon$, and let $\|g_i-\tilde g_i \|_{L^2(\Omega)} < \delta_i$ $(i=1,\dots ,r)$. Denote by $u^\ast = (u^\ast_1,\dots , u^\ast_r)$ and $\tilde{u} = (\tilde{u}_1,\dots , \tilde{u}_r)$ the solution s of the systems $$ T_i(u)= g_i(x) \quad \hbox{ and } \quad T_i(u)= \tilde g_i(x) $$ (both with boundary condition ${Qu_i}\big|_{\partial \Omega }=0$), respectively. Then $\|\tilde u-u^\ast\|_{H^1_Q(\Omega)^r} < \varepsilon$. \end{proposition} \begin{proof} Note that \begin{align*} m\|u^\ast-\tilde{u}\|_{H^1_Q(\Omega)^r}^2 &\le \sum_{i=1}^r \int_\Omega \bigl( T_i(u^\ast)-T_i(\tilde{u}) \bigr) (u^\ast_i-\tilde{u}_i )\\ & \le \sum_{i=1}^r \| g_i - \tilde g_i \|_{L^2(\Omega)} \| u^\ast_i-\tilde{u}_i \|_{L^2(\Omega)}\\ & \le \lambda_1^{-1/2} \sum_{i=1}^r \delta_i \| u^\ast_i-\tilde{u}_i \|_{H^1_Q} \\ &\le \lambda_1^{-1/2} (\sum_{i=1}^r \delta_i^2)^{1/2} \| u^\ast_i-\tilde{u}_i \|_{H^1_Q(\Omega)^r} < m\varepsilon \| u^\ast_i-\tilde{u}_i \|_{H^1_Q(\Omega)^r}\, . \end{align*} \end{proof} In the sequel we consider (\ref{ujSY}) as having polynomial coefficients, i.e. being replaced by the approximate system. For simplicity of presentation, we only consider the two-dimension al case (the analogies being straightforward). \subsection*{(a) Rectangle} We investigate the case when $\Omega\subset\mathbb{R}^2$ is a rectangle. It can be assumed that $\Omega=I\equiv [0,\pi]^2$ (if not, a linear transformation is done). Denote by $\mathcal{P}_s$ and $\mathcal{P}_c$ the set of sine- and cosine-polynomials \begin{gather*} \mathcal{P}_s= \{ \sum_{n,m=1}^l \sigma_{nm} \sin nx \sin my: \, l\in\mathbb{N}^+, \, \sigma_{nm}\in\mathbb{R} \, (n,m=1,\dots ,l) \},\\ \mathcal{P}_c= \{ \sum_{n,m=0}^l \sigma_{nm} \cos nx \cos my: \, l\in\mathbb{N}, \, \sigma_{nm}\in\mathbb{R} \, (n,m=0,\dots ,l) \}, \end{gather*} respectively. The coefficients $a_i$ and $c_j^i$ of (\ref{ujSY}) can be approximated by cosine-polynomials to any prescribed accuracy, hence (as suggested above) we assume that (\ref{ujSY}) already fulfils $a_i, c_j^i \in \mathcal{P}_c$ $(i=1,\dots ,r; |j|\le s_i)$. {\it Dirichlet boundary conditions.} We consider the case when $|j|$ is odd in each term in (\ref{ujSY}) and assume $g_i\in \mathcal{P}_s$. Then the operators $T_i$ are invariant on $\mathcal{P}_s$. Further, for any $h\in \mathcal{P}_s$ the solution of $$ -\Delta z = h, \quad z\big|_{\partial I}= 0 $$ fulfils $z\in \mathcal{P}_s$, namely, if $ h(x,y)= \sum_{n,m=1}^l \sigma_{nm} \sin nx \sin my$, then \begin{equation} z(x,y)= \sum_{n,m=1}^l {\sigma_{nm} \over n^2+m^2} \sin nx \sin my\, . \label{laplinv} \end{equation} These imply that if (\ref{iter-a}) starts with $u^0_i\in \mathcal{P}_s$ then we have $u^k_i\in \mathcal{P}_s$ throughout the iteration, and the auxiliary problems are solved in the trivial way (\ref{laplinv}). {\it Neumann boundary conditions.} Similar to the Dirichlet case, now letting $g_i\in\mathcal{P}_c$. Here $T_i$ are invariant on $\mathcal{P}_c$. Further, for $h\in \mathcal{P}_c$, the solution $ h(x,y)= \sum_{n,m=0}^l \sigma_{nm} \cos nx \cos my$, of $$ -\Delta z + cz = h, \quad \partial_\nu z\big|_{\partial I}= 0 $$ fulfils $z\in \mathcal{P}_c$, namely, $ z(x,y)= \sum_{n,m=0}^l {\sigma_{nm} \over n^2+m^2+c}\cos nx \cos my$. {\it Mixed boundary conditions.} Denote by $\Gamma_i$ $(i=1, \dots ,4)$ the boundary portions $[0,\pi)\times \{ 0 \}$, $\{ \pi \}\times [0,\pi)$, $(0,\pi]\times \{ \pi \}$, $\{ 0 \}\times (0,\pi]$, respectively. First, let $\alpha (x)\equiv \chi_{ \{ \Gamma_2 \cup \Gamma_4 \} }$, $\beta (x)\equiv \chi_{ \{ \Gamma_1 \cup \Gamma_3 \} }$, i.e. we have $ u\big|_{\Gamma_2 \cup \Gamma_4}\equiv 0$, $ \partial_\nu u\big|_{\Gamma_1 \cup \Gamma_3} \equiv 0$. Then the above method works on $ \mathcal{P}_{sc} = \{ \sum_{n=1}^l \sum_{m=0}^l \sigma_{nm} \sin nx \cos my \, \}$. We can proceed similarly for other edgewise complementary characteristic functions $\alpha$ and $\beta$, using $\sin (n+{1 \over 2})x$ type terms for mixed endpoint conditions. \subsection*{(b) Disc} We investigate the case when $\Omega\subset\mathbb{R}^2$ is the unit disc $B=B_1({\bf 0})$. Now $a_i$, $c_j^i$ and $g_i$ are assumed to be algebraic polynomials. (Notation: $a_i, c_j^i, g_i\in\mathcal{P}_{alg}$). Then $T_i$ are invariant on $\mathcal{P}_{alg}$. {\it Dirichlet boundary conditions.} If $h\in\mathcal{P}_{alg}$, then the solution of %the auxiliary problem \ $$ -\Delta z = h, \quad z\big|_{\partial B}= 0 $$ can be found by looking for $z$ in the form $$ z(x,y)=(x^2+y^2-1)q(x,y), $$ where $q\in\mathcal{P}_{alg}$ has the same degree as $h$ (cf. \cite{K-A}). Then the coefficients of $q$ can be determined by solving a a simply structured linear system of algebraic equations (obtained from equating the coefficients of $-\Delta z$ and $h$). The matrix of this system has three diagonals and two other non-zero elements in each row. \noindent{\it $3^{rd}$ boundary conditions.} We examine the case $\alpha(x)\equiv \alpha >0$, $\beta(x)\equiv \beta >0$. If $h\in\mathcal{P}_{alg}$ then the solution of $$ -\Delta z = h, \quad (\alpha z + \beta \partial_\nu z)\big|_{\partial B}= 0 $$ can be determined similarly to the Dirichlet case. Namely, let $q(x,y)$ be a polynomial with unknown coefficients $\{ a_{nm} \}$ and of the same degree as $h$. Let $$ p(x,y)\equiv (x^2+y^2-1)q(x,y) = \sum_{s=0}^l \sum_{n+m=s} b_{nm} x^n y^m, $$ here the $b_{nm}$'s are linear combination s of the $a_{nm}$'s. We look for $z$ as $$ z(x,y)= \sum_{s=0}^l \sum_{n+m=s} {b_{nm} \over {\alpha+\beta s}} x^n y^m\, . $$ Equating the coefficients of $-\Delta z$ and $h$ leads again to a linear system of algebraic equations for $\{ a_{nm} \}$, from which we then determine $\{ b_{nm} \}$. Then $z$ fulfils the boundary condition, since on $\partial B$ we have $\partial_\nu z = x \partial_x z + y \partial_y z$ and $$ \alpha z + \beta (x \partial_x z + y \partial_y z) = \sum_{s=0}^l \sum_{n+m=s} {b_{nm} \over {\alpha+\beta s}}(\alpha + \beta (n+m)) x^n y^m = p(x,y) = 0\, . $$ \subsection*{(c) Annulus} Let $A= \{ (x,y)\in\mathbb{R}^2: \, R_1^2 < x^2+y^2 < R_2^2\, \}$ with given $R_2>R_1>0$. Let $\Gamma_m = C({\bf 0}, R_m)$ $(m=1,2)$ be the boundary circles. We investigate the system with radially symmetric coefficients, written in the form \begin{equation} \begin{gathered} T_i(u)\equiv -\overline{\nabla} \bigl( a_i(x^2+y^2) \overline{\nabla} u_i \bigr) + \sum_{|j|\le s_i} c^{(i)}_j (x^2+y^2) u_1^{j_1}\dots u_n^{j_n} = g_i(x^2+y^2) \\ Qu_i\equiv (\alpha u_i +\beta \partial_\nu u_i) \big|_{\Gamma_m }=0 \quad (m=1,2) \end{gathered}\label{rad} \end{equation} ($i=1, \dots ,n$), with the notation $\overline{\nabla} u = \overline{\partial}_x u + \overline{\partial}_y u$, $\overline{\partial}_x u = x\partial_x u$, $\overline{\partial}_y u = y\partial_y u$. The functions $a_i\in C^1[R_1,R_2]$, $c^{(i)}_j\in C[R_1,R_2]$, $g_i\in L^2[R_1,R_2]$ and the numbers $\alpha, \beta\ge 0$ are such that the positivity and monotonicity conditions (C2)-(C4) are fulfilled. Introducing the notations $$ r=(x^2+y^2)^{1 \over 2}, \quad \hat a_i (r)= r a_i(r^2), \quad \hat c^{(i)}_j (r) = {1 \over r} c^{(i)}_j(r^2), \quad \hat g_i (r) = {1 \over r} g_i (r^2) $$ and using $\overline{\nabla} u = r \partial_r u$, the system (\ref{rad}) is written as \begin{align*} \hat T_i(u)\equiv - \frac{1}{2\pi r} \Big( \partial_r ( \hat a_i(r) \partial_r u_i ) + \sum_{|j|\le s_i} \hat c^{(i)}_j (r) u_1^{j_1}\dots u_n^{j_n} \Big) = \frac{1}{2\pi r} \hat g_i(r) \\ \hat Q_m u_i\equiv (\alpha u_i + (-1)^m \beta \partial_r u_i) \big|_{r=R_m}=0 \quad (m=1,2). \end{align*} Let $$ Bu\equiv \frac{1}{2\pi r} ( -\partial_r^2 u + cu ) \quad \hbox{on } D(B)=D(T_i)= H^2_{rad}(A) $$ where $H^2_{rad}(A)$ = \{$u\in H^2(A): u$ is radially symmetric and $Qu\big|_{\partial A}=0$ in trace sense\}, and $c=\kappa /m$ as in Section 2. Then there holds $$ \int_A (Bu)v \, dx\,dy = \int_{R_1}^{R_2} \Big( \partial_r u\, \partial_r v \, + cuv \Big)\, dr + \beta_Q\, \quad (u,v\in H^2_{rad}(A)), $$ where $\beta_Q =0$ if $\beta =0$ and $\beta_Q = {\alpha \over \beta}\sum_{i=1}^2 u(R_m)v(R_m)$ if $\beta>0$. Hence \, $H_B=H^1_Q(A)\subset$\{ $u\in H^1_0(A)$: $u$ is radially symmetric \}. Further, we have $$ \int_A \hat T_i (u)v \, dxdy = \int_{R_1}^{R_2} \Big( \hat a_i \partial_r u\partial_r v + \sum_{|j|\le s_i} \hat c^{(i)}_j u_1^{j_1}\dots u_n^{j_n}\, v \Big)\, dr + \tau^{(i)}_Q , $$ where $\tau^{(i)}_Q =0$ if $\beta =0$ and $\tau^{(i)}_Q = {\alpha \over \beta} \sum_{i=1}^2 \hat a_i(R_m)u(R_m)v(R_m)$ if $\beta>0$. Consequently, Theorem \ref{th21} applies in the setting $N=1$ on $(R_1,R_2)$. That is, (\ref{rad}) has a unique weak solution which is radial: $u^\ast=(u_1^\ast, \dots , u_n^\ast)\in H^1_Q(A)^n$. Further, for any $u_0\in H^2_{rad}(A)^n$ the sequence (\ref{iter-c}) converges linearly to $u^\ast$. The auxiliary problem (\ref{iter-c}) is equivalent to $$ -\partial_r^2 z_i^k = \hat g^k_i\, , \quad (\hat Q_m u_i) \big|_{r=R_m}=0 $$ where $\hat g^k_i \equiv -\partial_r (\hat a_i \partial_r u^k_i) + \sum_{|j|\le s_i} \hat c^{(i)}_j (u^k_1)^{j_1}\dots (u^k_n)^{j_n} - \hat g_i$. If $\alpha\not\equiv 0$ then $\hat a_i$, $\hat c^{(i)}_j$ and $\hat g_i$ are replaced by approximating algebraic polynomials. Then $\hat T_i$ are invariant on $\mathcal{P}_{rad}\equiv$\{ algebraic polynomials of $r$\}. For any $h\in\mathcal{P}_{rad}$ the solution of $$ -\partial_r^2 z = h, \quad (\hat Q_m z) \big|_{r=R_m}=0 $$ also fulfils $z\in\mathcal{P}_{rad}$ and is elementary to determine. Namely, if $h(r)=\sum_{j=0}^l a_j (r-R_1)^j$, then $$ z(r)=-\sum_{j=0}^l \frac{a_j}{(j+1)(j+2)} (r-R_1)^{j+2} \, + c(r-R_1) + d $$ where the constants $c$ and $d$ are determined from the boundary condition. In the case $\alpha\equiv 0$ the iteration is kept in the class of \ cosine-polynomials of $r$ on $(R_1,R_2)$ and the auxiliary problems are solved as in subsection (a) with Neumann conditions. \subsection*{(d) Other domains} The methods given in paragraph (a) for a rectangle can be extended for other domains where the eigenfunctions of the Laplacian are known explicitly. Then the terms of the polynomials in $\mathcal{P}_s$, $\mathcal{P}_c$ and $\mathcal{P}_{sc}$ have to be replaced by the actual eigenfunctions. These are known e.g. for rectangular triangles and regular hexagons \cite{Makai, R-Sz}. \subsection*{(e) Domain transformation} If two domains are diffeomorphic, then (\ref{ujSY}) can be transformed from one to the other such that uniform monotonicity is preserved. The formulation of this property is restricted to the Dirichlet problem. Further, for simplicity of notation it is done for a single equation. \begin{proposition} \label{prop52} Let $S$, $\Omega\subset\mathbb{R}^N$ be bounded domains, $\Phi\in C^1(\overline S, \overline\Omega)$, $\det \Phi'(x)\neq 0$ $(x\in \overline S)$. Let $F:H^1_0(\Omega)\to H^1_0(\Omega)$ be given by $$ \langle F(u),v\rangle _{H^1_0(\Omega)} = \int_\Omega \bigl( a(x) \, \nabla u \cdot \nabla v + \sum_{j=0}^s c_j(x)u^j v \bigr)\, dx \quad (u, v\in H^1_0(\Omega)) $$ where $a$ and $c_j$ fulfil conditions (C1)-(C4), which means here that $a\in C^1(\overline\Omega)$, $ c_j\in C(\overline\Omega)$, $a(x)\ge m>0$, $0\le \sum_{j=0}^s j c_j(x)u^{j-1}$, $s\le p-1$ with $p$ in (C4). (That is, $F$ is the generalized operator corresponding to some differential operator $T$ of the studied kind.) Let $\tilde{u}\equiv u\circ \Phi$ $(u\in L^2(\Omega))$. Then the operator $\tilde F:H^1_0(S)\to H^1_0(S)$, defined by $$ \langle \tilde F(\tilde u),\tilde v\rangle _{H^1_0(S)}= \langle F(u),v\rangle _{H^1_0(\Omega)} \quad (u,v\in H^1_0(\Omega)), $$ fulfils $$ m_1 \|\tilde h \|_{H^1_0(S)}^2 \le \langle \tilde F'(\tilde u)\tilde h,\tilde h\rangle _{H^1_0(S)} \le\, M_1(\tilde{u}) \|\tilde h\|_{H^1_0(S)}^2 \quad (\tilde u,\tilde h\in H^1_0(S)) $$ with suitable constants $m_1$ and $M_1(\tilde{u})$, the latter depending on $\tilde{u}$. \end{proposition} \begin{proof} Setting $\tilde a = (a\circ \Phi)\det \Phi'$ and $\tilde c_j = (c_j\circ \Phi){\rm det}\Phi'$, we have $$ \langle \tilde F(\tilde u),\tilde v\rangle _{H^1_0(S)}= \int_S \bigl( \tilde a \, \widetilde{\nabla u} \cdot \widetilde{\nabla v} + \sum_{j=0}^s \tilde c_j \tilde u^j \tilde v \bigr)\, dx \, . $$ Then there holds $$ \langle \tilde F'(\tilde u)\tilde h,\tilde h\rangle _{H^1_0(S)}= \int_S \bigl( \tilde a |\widetilde{\nabla h}|^2 + \sum_{j=0}^s j \tilde c_j \tilde u^{j-1} \tilde h^2 \bigr)\, dx \, . $$ Using $0< \min_{\overline\Omega} {\rm det}\Phi' $ and $\max_{\overline\Omega} {\rm det}\Phi'< +\infty$, it is easy to verify that $\tilde F'(\tilde u)$ inherits uniform ellipticity from $F'(u)$. The details of calculations are left to the reader. \end{proof} If $\Phi\in C^2$, then, due to the smoothness conditions and the ellipticity result of Proposition \ref{prop52}, the differential operator $\tilde T$ that the generalized operator $\tilde F$ corresponds to inherits the properties of the original operator $T$. Consequently, if $\Omega$ is $C^2$-diffeomorphic to one of the above special domains (say $S$), then transforming the equation $T(u)=g$ from $\Omega$ to $S$ via $\Phi$ leads to the described direct realization. Finally, we note that the method described in this section also works for the general form of our system (\ref{SY}) if the nonlinearity $f(x,u)$ can be suitably approximated by polynomials, leading to the form (\ref{ujSY}). \section{Examples} % 3.3 \label{secex} We consider two model problems on the domain $I=[0,\pi]^2\subset \mathbb{R}^2$: (1) a single equation with Dirichlet boundary conditions; (2) a system of two equations with mixed boundary conditions. According to section \ref{secdirect} (a), the operations required for the numerical realization are elementary. Namely, the terms of the iteration are trigonometric polynomials. The iteration simply consists of executing multiplications and additions of the polynomials, further (for $-\Delta^{-1}$), linear combination of the form (\ref{laplinv}). Storing the polynomials as matrices of the coefficients, the algorithmization is straightforward. \subsection*{General considerations and the algorithm} There is one more step to be inserted in the realization. Namely, the exact execution of the iteration yields exponentially growing degree of the polynomials. Trying to avoid this for obvious memory reasons, it turns out that the amount of terms which would cause the storing problem consists of essentially useless (almost 0) terms (the high-index coefficients of the stored polynomials), and most of them can be neglected with small prescribed change of accuracy. We formulate this below generally and conclude at the final algorithm for our example using a simple truncation procedure. (The idea may be put through to the other special domains.) \begin{definition} \label{def61} \rm Let $\lambda_{kl}$ and $u_{kl}$ $(k,l\in\mathbb{N}^+)$ be the eigenvalues and eigenfunctions (normed in $L^2(\Omega)$) of the auxiliary problem $$ -\Delta u + cu = \lambda u, \quad Qu\big|_{\partial \Omega} = 0\, . $$ Let $u\in H^1_Q(\Omega)$ be fixed, $u=\sum_{k,l=1}^{\infty} a_{kl} u_{kl}$. Then \begin{enumerate} \item[(a)] for any $s\in\mathbb{N}^+$ we define $ u_{[s]}\equiv \sum_{k+l\le s} a_{kl} u_{kl}$; \item[(b)] for any $0<\omega < \|u\|_{H^1_Q}$ the index $k_{u,\omega}\in\mathbb{N}^+$ is defined by the inequalities \begin{equation} \sum_{k+l> k_{u,\omega} } \lambda_{kl} a_{kl}^2 \le \omega^2, \quad \sum_{k+l\ge k_{u,\omega} } \lambda_{kl} a_{kl}^2 > \omega^2 .\label{KD} \end{equation} Since $\|u\|_{H^1_Q}^2=\sum_{k+l\ge 1} \lambda_{kl} a_{kl}^2$, the index $k_{u,\omega}\in\mathbb{N}^+$ is the smallest one for which \begin{equation} \|u-u_{[k_{u,\omega}]} \|_{H^1_Q}\le \omega\, . \label{szeletbecs} \end{equation} If the series of $u$ consists of finitely many terms itself, then we can also define $k_{u,\omega}$ via (\ref{KD}) for $\omega=0$. \end{enumerate} \end{definition} \begin{remark} \label{rem61} \rm In the considered examples we have $\lambda_{kl}= k^2+l^2$. \end{remark} \begin{proposition} \label{prop61} Let $0<\omega < \|u\|_{H^1_Q}$ be fixed, $u\in H^1_Q(\Omega)$, $(u^n)_{n\in\mathbb{N}}\subset H^1_Q(\Omega)$. If $u^n\to u$ in $H^1_Q(\Omega)$ then the sequence $(k_{u^n, \omega})$ is bounded. \end{proposition} \begin{proof} Let $n_0\in\mathbb{N}^+$ such that $\|u^n-u\|_{H^1_Q(\Omega)} \le {\omega \over 2}$ $(n\ge n_0)$. Let $u=\sum_{k,l=1}^{\infty} a_{kl} u_{kl}$, $u^n =\sum_{k,l=1}^{\infty} a_{kl}^n u_{kl}$, and $k_2= k_{u,\omega /2}$. Then for $n\ge n_0$ we have \begin{align*} \sum_{k+l> k_2} \lambda_{kl} (a^n_{kl})^2 &\le 2 \Big(\sum_{k+l> k_2} \lambda_{kl} (a^n_{kl}-a_{kl})^2 + \sum_{k+l> k_2} \lambda_{kl} a_{kl}^2 \Big) \\ &\le 2 \left((\omega / 2)^2 + (\omega / 2 )^2 \right) = \omega^2 . \end{align*} Hence $k_{u^n, \omega} \le k_2$ $(n\ge n_0)$. \end{proof} Note that the property holds similarly in product space. This proposition means that we can consider bounded number of terms in the series of $u^n$ to attain $u$ within a prescribed error. This motivates the following truncation procedure in the realization of the gradient method for our model problems on $I$. {\it The algorithm.} Let the coefficients of $T_i$ and the right sides $g_i$ be trigonometric polynomials as in section \ref{secdirect} (a) (after suitable approximation). The sequence $u^n=(u^n_1, \dots , u^n_r)\in H^1_Q(I)^r$ is constructed as follows. Let $u^0=(0, \dots , 0)$. If, for $n\in\mathbb{N}$, $u^n$ is obtained, then for $i=1,\dots ,r$ we define \begin{equation} \label{ALG} \begin{aligned} & g_i^n= T_i(u^n)-g_i, \\ & z_i^n \quad \hbox {is the solution of} \quad -\Delta z_i^n = g_i^n, \quad Qz_i^n \big|_{\partial I}=0; \\ &\hbox{ we fix } \omega_n>0,\\ &\tilde z_i^n = (z_i^n)_{[k_{z_i^n , \omega_n}]} \, ,\\ &u_i^{n+1} = u_i^n - \frac{2}{M+m} \tilde z_i^n\, . \end{aligned} \end{equation} \begin{proposition} \label{prop62} Let $\varepsilon >0$, $\omega=r^{-1/2}m\varepsilon$ (where $m$ is the lower bound of $T$), $(\omega_n)\subset\mathbb{R}^+$, $\omega_n=\omega$ $(n\ge n_0)$. Then there exists $c>0$ such that \ \begin{equation} \label{p36} \|u^n-u^\ast \|_{H^1_Q(I)^r} \le \varepsilon + c\big( \frac{M-m}{M+m} \big)^n \quad (n\in\mathbb{N}). \end{equation} \end{proposition} \begin{proof} Estimate (\ref{szeletbecs}) implies $\|z^n-\tilde z^n \|_{H^1_Q(I)^r} \le r^{1/2}\omega$ $(n\ge n_0)$, hence we can apply Corollary \ref{cor42} with initial guess $u^{n_0}$. \end{proof} In each step we calculate the residual errors $$ r_i^n = \|T_i(u^n)-g_i\|_{L^2(I)}\, . $$ Since $T_i(u^n)$ and $g_i$ are trigonometric polynomials, this only requires square summation of the coefficients. Letting $$ e_i^n = \frac{1}{m\sqrt{\lambda}} r_i^n\, , $$ we then obtain \begin{equation} e^n \equiv \, \|u^n-u^\ast\|_{H^1_Q(I)^r} \le \Big( \sum_{i=1}^r (e_i^n)^2 \Big)^{1/2} . \label{EN} \end{equation} \begin{remark} \label{rem63} \rm (a)\; Since $z_i^n$ is a trigonometric polynomial, therefore (writing $z_i^n = \sum_{k+l\le s_n} \zeta_{kl} a_{kl}$ and $k_n = k_{z^n_i, \omega_n}$) \, (\ref{KD}) means a condition on finite sums: $$ \sum_{k_n < k+l \le s_n} (k^2+l^2) \zeta_{kl}^2 \le \omega^2 , \quad \sum_{k_n \le k+l \le s_n} (k^2+l^2) \zeta_{kl}^2 > \omega^2 , $$ i.e. $\tilde z_i^n$ is determined by calculating appropriate square sums of cross-diagonals in the coefficient matrices. \noindent(b)\; The choice of the constants $\omega_n$ is influenced by the following factors. The eventual values of $\omega_n$ (i.e. for large $n$) are determined by the required accuracy in virtue of Proposition \ref{prop62}. To keep the matrix sizes low, we choose $\omega_n$ large as long as possible. Based on (\ref{p36}), the necessity for decreasing $\omega_n$ arises when the error $e^n$ (calculated by (\ref{EN})) ceases to follow the theoretical linear convergence $( \frac{M-m}{M+m})^n$. \end{remark} The truncation of the polynomials can be analoguously inserted in the algorithm in the case of other special domains. \begin{example} \rm % Example 1 We consider the single Dirichlet problem \begin{gather*} -\Delta u + u^3 = \sin x \sin y \\ u\big|_{\partial I }=0. \end{gather*} Since $g(x,y)=\sin x \sin y \in \mathcal{P}_s$, therefore, defining $u^0=0$, we have $(u^n)\subset\mathcal{P}_s$. According to the previous sections, the realization of the gradient method means elementary matrix operations for the multiplication and addition of sine polynomials (stored as matrices of coefficients), further, linear combination of the form (\ref{laplinv}) for solving the auxiliary problems \begin{gather*} -\Delta z^n = g^n \\ {z^n}\big|_{\partial I }=0 \end{gather*} (with $g^n= -\Delta u^n + (u^n)^3 -g$). The calculations are made up to accuracy $10^{-4}$. First the constants $m$ and $M$ for (\ref{ALG}) are determined from (\ref{Mdef}) in Theorem \ref{th21}. We now have $m=m'=1$, $\eta=0$, $p=4$, $\kappa=\kappa'=0$, $\gamma =3$, $\lambda_1=1^2+1^2=2$, $K_{4,I}\le 1$ from Corollary \ref{cor25} (a), $\mu_4=1$, and $\|g\|_{L^2(I)}={\pi \over 2}$. Hence we have $$ M=4.7011 \quad \hbox{ and thus } \quad \frac{2}{M+m} = 0.3508. $$ The theoretical convergence quotient is $$ \frac{M-m}{M+m} = 0.6492\, . $$ The algorithm (\ref{ALG}) has been performed in MATLAB, which is convenient for the matrix operations. The following table contains the error $e^n$ (see (\ref{EN})) versus the number of steps $n$. The truncation constant is $\omega_n \equiv \omega= 0.0005$ throughout the iteration. \begin{table}[ht] \noindent\begin{tabular}{|c|c|c|c|c|c|c|c|c|}\hline step $n$ & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\ \hline {error $e^n$} & 1.1107 & 0.7186 & 0.4555 & 0.2821 & 0.1717 & 0.1034 & 0.0620 & 0.0375 \\ \hline \end{tabular} \\[3pt] \begin{tabular}{|c|c|c|c|c|c|c|c|c|}\hline step $n$ & 9 & 10 & 11 & 12 & 13 & 14 & 15 & 16 \\ \hline error $e^n$ & 0.0231 & 0.0149 & 0.0093 & 0.0058 & 0.0036 & 0.0022 & 0.0014 & 0.0009 \\ \hline \end{tabular} \\[3pt] \begin{tabular}{|c|c|c|c|c|c|c|}\hline step $n$ & 17 & 18 & 19 & 20 & 21 & 22 \\ \hline error $e^n$ & 0.0006 & 0.0004 & 0.0003 & 0.0002 & 0.0002 & 0.0001 \\ \hline \end{tabular} %\caption{} %1 \end{table} \end{example} \begin{example} \rm %Example 2 We consider the system with mixed boundary conditions \begin{equation} \label{EX2} \begin{gathered} -\Delta u + u - v + u^3 = g_1(x,y) \\ -\Delta v + v - u + v^3 = \quad 0 \phantom{xy} \\ u\big|_{\Gamma_1}= v\big|_{\Gamma_1}=0, \partial_\nu u\big|_{\Gamma_2}= \partial_\nu v\big|_{\Gamma_2}=0 \end{gathered} \end{equation} where $$ g_1(x,y)= \frac{\sin x \cos y}{(2-0.249\cos 2x)(2-0.249\cos 2y)} $$ and $\Gamma_1=\{ 0,\pi \}\times [0,\pi]$, $\Gamma_2= [0,\pi ]\times \{ 0,\pi \}$. Now the eigenfunctions are $\sin kx \cos ly$ which are in $\mathcal{P}_{sc}$; hence $g_1$ is approximated by $\tilde g_1\in \mathcal{P}_{sc}$. Calculating up to accuracy $10^{-4}$, we define the Fourier partial sum $$ \tilde g_1 (x,y)= \sum_{k,l are odd \atop k+l\le 6} a_{kl} \sin kx \cos ly\, , \quad a_{kl}=4.0469\cdot 4^{-(k+l)} $$ which fulfils $\|g_1-\tilde g_1 \|_{L^2(I)} \le 0.0001$. Replacing $g_1$ by $\tilde g_1$ in (\ref{EX2}), Proposition \ref{prop51} yields $$ \|\tilde u-u^\ast\|_{H^1_Q(I) ^2} \le 0.0001, $$ where $u^\ast = (u^\ast_1,u^\ast_2)$ and $\tilde{u} = (\tilde{u}_1,\tilde{u}_2)$ are the solution s of the original and approximated systems, respectively. Defining $u^0=v^0=0$, we have $(u^n)\subset \mathcal{P}_{sc}$ and $(v^n)\subset \mathcal{P}_{sc}$. The iteration of the gradient method is realized through elementary matrix operations in an analogous way to Example 1. The calculations are made up to accuracy $10^{-4}$ again. The constants for (\ref{Mdef}) are determined as follows. We have $m=m'=1$, $\eta=0$, $p=4$, $\lambda_1=1^2+0^2=1$, $\mu_4=1$ from the given functions. Further, the eigenvalue s of the Jacobian of $f(u,v)=(u-v+u^3,v-u+v^3)$ are between 0 and $2+3(u^2+v^2)$, hence $\kappa'=2$ and $\gamma =3$. Besides, we can use Lemmas 2.1--2.2: the boundary condition gives $K_{2,\Gamma_1}=0$ and Lemma \ref{lem22} yields $K_{2,\Gamma_2}\le 2.1535$, hence by Lemma \ref{lem21} we have $K_{4,I}\le 1.6051$. Finally, $\|\tilde g_1\|_{L^2(I)}=0.3988$. Hence we have $$ M=6.1420 \quad \hbox{and thus} \quad \frac{2}{M+m} = 0.2801, $$ the theoretical convergence quotient now being $$ \frac{M-m}{M+m} = 0.7200\, . $$ The following table presents the errors (first given for $u^n$ and $v^n$, resp.). The truncation constant is $\omega_n \equiv \omega= 0.0005$ up to step 10, then (observing that the error quotient rises over 0.8) $\omega_n \equiv \omega= 0.0001$ is chosen. \begin{table}[ht] \begin{tabular}{|c|c|c|c|c|c|c|c|}\hline step $n$ & 1 & 2 & 3 & 4 & 5 & 6 & 7 \\ \hline $e^n_1$ & 0.2821 & 0.1219 & 0.0627 & 0.0358 & 0.0215 & 0.0132 & 0.0083 \\ \hline $e^n_2$ & 0 & 0.0532 & 0.0459 & 0.0316 & 0.0204 & 0.0129 & 0.0081 \\ \hline $[(e^n_1)^2 + (e^n_2)^2 ]^{1/2}$ & 0.2821 & 0.1330 & 0.0776 & 0.0477 & 0.0297 & 0.0184 & 0.0116 \\ \hline \end{tabular} \\[3pt] \begin{tabular}{|c|c|c|c|c|c|c|c|}\hline step $n$ & 8 & 9 & 10 & 11 & 12 & 13 & 14 \\ \hline $e^n_1$ & 0.0053 & 0.0036 & 0.0027 & 0.0022 & 0.0014 & 0.0009 & 0.0005 \\ \hline $e^n_2$ & 0.0053 & 0.0038 & 0.0030 & 0.0025 & 0.0015 & 0.0010 & 0.0006 \\ \hline $[ (e^n_1)^2+ (e^n_2)^2 ]^{1/2}$ & 0.0075 & 0.0052 & 0.0041 & 0.0033 & 0.0020 & 0.0013 & 0.0008 \\ \hline \end{tabular} \\[3pt] \begin{tabular}{|c|c|c|c|c|c|c|c|}\hline step $n$ & 15 & 16 & 17 & 18 \\ \hline $e^n_1$ & 0.0003 & 0.0002 & 0.0001 & 0.0001 \\ \hline $e^n_2$ & 0.0003 & 0.0002 & 0.0001 & 0.0001 \\ \hline $[ (e^n_1)^2 + (e^n_2)^2 ]^{1/2}$ & 0.0005 & 0.0003 & 0.0002 & 0.0001 \\ \hline \end{tabular} %\caption{} %2 \end{table} \end{example} \begin{thebibliography}{00} \bibitem{Ad} {\sc Adams, R.A.}, {\it Sobolev spaces}, Academic Press, New York-London, 1975. \bibitem{Ax} {\sc Axelsson, O.,} {\it Iterative solution methods,} Cambridge Univ. Press, 1994. \bibitem{Ax-C} {\sc Axelsson, O., Chronopoulos, A.T.,} {\it On nonlinear generalized conjugate gradient methods,} Numer. Math. 69 (1994), No. 1, pp. 1-15. \bibitem{Ax-G} {\sc Axelsson, O., Gustafsson, I.,} {\it An efficient finite element method for nonlinear diffusion problems,} Bull. GMS, 32 (1991), pp. 45-61. \bibitem{Ax-Lay} {\sc Axelsson, O., Layton, W.,} {\it Iteration methods as discretiz ation procedures,} Lecture Notes in Math. 1457, Springer, 1990. \bibitem{BW} {\sc B\"orgers, C., Widlund, O. B.,} On finite element domain imbedding methods, {\it SIAM J. Numer. Anal.} 27 (1990), no. 4, 963--978. \bibitem{Bur1} {\sc Burenkov, V.I., Gusakov, V.A.}, {\it On exact constants in Sobolev embeddings} (in Russian), Dokl. Akad. Nauk. 294 (1987), No. 6., pp. 1293-1297. \bibitem{Bur2} {\sc Burenkov, V.I., Gusakov, V.A.}, {\it On exact constants in Sobolev embeddings III.}, Proc. Stekl. Inst. Math. 204 (1993), No. 3., pp. 57-67. \bibitem{Da} {\sc Daniel, J.W.}, {\it The conjugate gradient method for linear and nonlinear operator equations}, SIAM J. Numer. Anal., 4 (1967), No.1., pp. 10-26. \bibitem{Dorr} {\sc Dorr, F. W.,} The direct solution of the discrete Poisson equation on a rectangle, {\it SIAM Rev.} 12 (1970), 248--263. \bibitem{Enc} {\sc Egorov, Yu.V., Shubin, M.A.}, {\it Encyclopedia of Mathematical Sciences, Partial Differential Equations I.}, Springer, 1992. \bibitem{book} {\sc Farag\'o I., Kar\'atson J.,} {\it Numerical solution of nonlinear elliptic problems via preconditioning operators: theory and application.} Advances in Computation, Volume 11, {\it NOVA Science Publishers}, New York, 2002. \bibitem{GFEM} {\sc Farag\'o, I., Kar\'atson, J.,} The gradient-finite element method for elliptic problems, {\it Comp. Math. Appl.} {\bf 42} (2001), 1043-1053. \bibitem{Funaro92} {\sc Funaro, D.,} {\it Polynomial approximation of differential equations,} Lecture Notes in Physics, New Series, Monographs 8, Springer, 1992. \bibitem{GGZ} {\sc Gajewski, H., Gr\"oger, K., Zacharias, K., }{\it Nichtlineare Operatorgleichungen und Operatordifferentialgleichungen,} Akademie-Verlag, Berlin, 1974 \bibitem{Hack} {\sc Hackbusch, W., }{\it Theorie und Numerik elliptischer Differentialgleichungen}, Teubner, Stuttgart, 1986. \bibitem{JAA} {\sc Kar\'atson, J., } The gradient method for non-differentiable operators in product Hilbert spaces and applications to elliptic systems of quasilinear differential equations, {\it J. Appl. Anal.}, { 3} (1997) No.2., pp. 205-217. \bibitem{cgmn} {\sc Kar\'atson, J., } The conjugate gradient method for a class of non-differentiable operators, {\it Annales Univ. Sci. ELTE}, {\bf 40} (1997), 121-130. \bibitem{nonloc} {\sc Kar\'atson J.,} Gradient method in Sobolev space for nonlocal boundary value problems, {\it Electron. J. Diff. Eqns.}, Vol. {\rm 2000} (2000), No. 51, pp. 1-17. \bibitem{varprec} {\sc Kar\'atson J., Farag\'o I.,} Variable preconditioning via quasi-Newton methods for nonlinear problems in Hilbert space, {\it SIAM J. Numer. Anal.} {\rm 41} (2003), No. 4, 1242-1262. \bibitem{elpot} {\sc Kar\'atson J., L\'oczi L.,} Sobolev gradient preconditioning for the electrostatic potential equation, to appear in {\it Comp. Math. Appl.} \bibitem{K-A} {\sc Kantorovich, L.V., Akilov, G.P., }{\it Functional Analysis, } Pergamon Press, 1982. \bibitem{Kelley} {\sc Kelley, C.T.,} {\it Iterative methods for linear and nonlinear equations,} Frontiers in Appl. Math., SIAM, Philadelphia, 1995. \bibitem{Li} {\sc Lions. P.-L., Pacella, F., Tricario, M.,}, {\it Best constants in Sobolev inequalities for functions vanishing on some part of the boundary}, Indiana Univ. Math. J., 37 (1988), No. 2., pp. 301-324. \bibitem{Makai} {\sc Makai, E.,} Complete systems of eigenfunctions of the wave equation in some special case, {\it Studia Sci. Math. Hung.,} 11 (1976), 139--144. \bibitem{O-R} {\sc Ortega, J.M., Rheinboldt, W.C.,} {\it Iterative solution of nonlinear equations in several variables,} Academic Press, 1970. \bibitem{Ne} {\sc Necas, J., }{\it Introduction to the Theory of Elliptic Equations, } J.Wiley and Sons, 1986. \bibitem{NeuS} {\sc Neuberger, J. W., Renka, R. J.}, {\it Sobolev gradients and the Ginzburg-Landau functional,} SIAM J. Sci. Comput. 20 (1999), No. 2, pp. 582--590. \bibitem{Neu1} {\sc Neuberger, J. W.}, {\it Sobolev gradients and differential equations}, Lecture Notes in Math., No. 1670, Springer, 1997. \bibitem{Neu3} {\sc Neuberger, J. W.}, {\it Steepest descent for general systems of linear differential equations in Hilbert space }, Lecture Notes in Math., No. 1032, Springer, 1983. \bibitem{R-Sz} {\sc Riesz F., Sz.-Nagy B.,} {\it Vorlesungen \"uber Funktionalanalysis}, Verlag H. Deutsch, 1982. \bibitem{Ros} {\sc Rosen, G.}, {\it Minimum value for $c$ in the Sobolev inequality $\Vert \phi\sp{3}\Vert \leq c\Vert \nabla \phi\Vert \sp{3}$}, SIAM J. Appl. Math., 21 (1971), pp. 30-32. \bibitem{Swarz77} {\sc Swarztrauber, P. N.,} The methods of cyclic reduction, Fourier analysis and the FACR algorithm for the discrete solution of Poisson's equation on a rectangle, {\it SIAM Rev.} 19 (1977), no. 3, 490--501. \bibitem{Tal} {\sc Talenti, G.}, {\it Best constants in Sobolev inequality}, Ann. Mat. Pura Appl., Ser. 4 (1976), Vol. 110, pp. 353-372. \end{thebibliography} \end{document}