---
abstract: 'In this paper we develop an $hp$-adaptive procedure for the numerical solution of general, semilinear elliptic boundary value problems in 1d, with possible singular perturbations. Our approach combines both a prediction-type adaptive Newton method and an $hp$-version adaptive finite element discretization (based on a robust [*a posteriori*]{} residual analysis), thereby leading to a fully $hp$-adaptive Newton-Galerkin scheme. Numerical experiments underline the robustness and reliability of the proposed approach for various examples.'
address:
- 'Lucerne University of Applied Sciences and Arts, CH-6002 Luzern, Switzerland'
- ' Institut für Analysis und Scientific Computing, TU Wien, A-1040 Wien, Austria'
- 'Mathematics Institute, University of Bern, CH-3012 Switzerland'
author:
- Mario Amrein
- 'Jens M. Melenk'
- 'Thomas P. Wihler'
bibliography:
- 'references.bib'
title: 'An $hp$-Adaptive Newton-Galerkin Finite Element Procedure for Semilinear Boundary Value Problems'
---

Introduction
============

The purpose of this paper is the development of a numerical approximation procedure for semilinear elliptic boundary value problems, with possible singular perturbations. More precisely, for a fixed parameter ${\epsilon}>0$ (possibly with ${\epsilon}\ll 1$), and a continuously differentiable function $f:\,\Omega\times\mathbb{R}\to\mathbb{R}$, we consider the problem of finding a function $u:\,\Omega\to\mathbb{R}$ in an interval $\Omega=(a,b)$, $a<b$, which satisfies $$\label{eq:poisson}
\begin{aligned}
-{\epsilon}u''(x) &=f(x,u(x)) \text{ for }x\in \Omega,\qquad
u(a)=u(b)=0. 
\end{aligned}$$ Solutions of  are typically not unique (even infinitely many solutions may exist), and, in the singularly perturbed case, may exhibit boundary layers, interior shocks, and (multiple) spikes. The existence of multiple solutions due to the nonlinearity of the problem and/or the appearance of singular effects constitute two challenging issues when solving problems of this type numerically; see, e.g.,[@RoStTo08; @Verhulst].

*Newton-Galerkin Methods:* In order to address the numerical solution of nonlinear singularly perturbed problems of type , a Newton-Galerkin approach has recently been proposed in [@AmreinWihler:14; @AmreinWihler:15]. It is built upon two main ideas: Firstly, adaptive Newton methods are applied on the continuous level in order to approximate the nonlinear differential equation  by a sequence of linearized problems. Secondly, the resulting linear problems are discretized by means of an adaptive finite element procedure which, in turn, is based on an ${\epsilon}$-robust *a posteriori* error analysis. Then, providing an appropriate interplay between the adaptive Newton method and the adaptive finite element approach, by which we either perform a Newton step (if the Newton linearization effect dominates) or refine the current finite element mesh based on some suitable [*a posteriori*]{} error indicators (in case that the finite element discretization causes the main source of error), a fully adaptive Newton-Galerkin scheme is obtained; see also the related works on monotone nonlinear problems [@ChaillouSuri:07; @El-AlaouiErnVohralik:11; @CongreveWihler:15], or the articles [@ChaillouSuri:06; @Han:94] dealing with modelling errors in linearized models. Our numerical results in [@AmreinWihler:15] reveal that sensible solutions can be found even in the singularly perturbed case, and that our scheme is reliable in the sense that the adaptive Newton-Galerkin method is able to converge to solutions which preserve the qualitative structure of (reasonably chosen) initial guesses (see also [@AmreinWihler:14; @ScWi11]).

*$hp$-FEM and $hp$-adaptivity:* The aim of the present paper is to extend our work in [@AmreinWihler:15] from a low-order to an $hp$-adaptive Newton-Galerkin approach. For simplicity of exposition, we focus on the 1d case, however, we mention that an extension to the multi-d case is possible as well. Exploiting the $hp$-framework may potentially lead to *exponential rates of convergence* in the numerical approximation. For the purpose of $hp$-adaptivity, we carry out an ${\epsilon}$-robust $hp$-version *a posteriori* residual analysis for the FEM discretizations of the Newton linearizations. We remark that using $hp$-adaptivity, in contrast to the adaptive low-order methodology developed in the previous works [@El-AlaouiErnVohralik:11; @AmreinWihler:14; @AmreinWihler:15], introduces an additional aspect of adaptivity that may further improve the performance of the Newton-Galerkin approach. In the $hp$-literature, several adaptive strategies and algorithms have been proposed; see, e.g., the overview article [@MiMc11] and the references therein. In particular, we point to the smoothness estimation techniques proposed in [@eibner-melenk06; @HoustonSeniorSuliENUMATH; @HoustonSuliHPADAPT; @Ma94], and the related approach [@FaWiWi14; @Wi11_2; @Wi11] involving continuous Sobolev embeddings. The latter works will be applied in the current paper in order to drive the $hp$-adaptive finite element refinements; in combination with the adaptive Newton linearizations this will lead to a fully $hp$-adaptive Newton-Galerkin procedure.

*Problem formulation and notation:* Henceforth, we suppose that a (not necessarily unique) solution $u\in X:=H^1_0(\Omega)$ of  exists; here, we denote by $H^1_0(\Omega)$ the standard Sobolev space of functions in $H^1(\Omega)=W^{1,2}(\Omega)$ with zero trace on $\partial\Omega:=\{a,b\}$. Furthermore, signifying by $X'=H^{-1}(\Omega)$ the dual space of $X$, and upon defining the map ${\mathsf{F}}_{{\epsilon}}: X\rightarrow X'$ through $$\label{eq:Fweak}
\left \langle {\mathsf{F}}_{{\epsilon}}(u),v\right \rangle :=  \int_{\Omega}\left\{{\epsilon}u'v'-f(u)v\right\}{\,{\mathsf{d}}x}\qquad \forall v\in X,$$ where $\left\langle\cdot,\cdot\right\rangle$ is the dual product in $X'\times X$, the above problem  can be written as a nonlinear operator equation in $X'$: $$\label{eq:F0}
u\in X:\qquad {\mathsf{F}}_{\epsilon}(u)=0.$$ In addition, for open subsets $D\subseteq\Omega$, we introduce the norm $${\left|\!\left|\!\left|u\right|\!\right|\!\right|}_{{\epsilon},D}:=\Bigl({\epsilon}{\left\|u'\right\|}_{0,D}^2 +{\left\|u\right\|}_{0,D}^2 \Bigr)^{\nicefrac{1}{2}},$$ where $\|\cdot\|_{0,D}$ denotes the $L^2$-norm on $D$. Note that, in the case of $f(u)=-u+g$, with given $g\in L^2(\Omega)$, i.e., when  is linear and strongly elliptic, the norm ${\left|\!\left|\!\left|\cdot\right|\!\right|\!\right|}_{{\epsilon},\Omega}$ is a natural energy norm on $X$. Frequently, for $D=\Omega$, the subindex ‘$D$’ will be omitted. Furthermore, the associated dual norm of ${\mathsf{F}}_{\epsilon}$ from  is given by $${\left|\!\left|\!\left|{\mathsf{F}}_{\epsilon}(u)\right|\!\right|\!\right|}_{X',{\epsilon}} =  \sup_{\genfrac{}{}{0pt}{}{v\in X}{{\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon}}=1}}\int_{\Omega}\left\{{\epsilon}u'v'-f(u)v\right\}{\,{\mathsf{d}}x}.$$ Throughout this work we shall use the abbreviation $x \preccurlyeq y$ to mean $x \leq c y$, for a constant $c>0$ independent of the mesh size $h$, the local polynomial degree $p$, and of ${\epsilon}>0 $.

*Outline:* In Section \[sc:anewton\] we will consider the Newton method within the context of dynamical systems in general Banach spaces, and present a prediction strategy for controlling the Newton step size parameter. Furthermore, Section \[sc:FEM\] focuses on the application of the Newton methodology to semilinear elliptic problems as well as on the discretization of the problems under consideration by $hp$-finite element methods; in addition, we derive an ${\epsilon}$-robust [*a posteriori*]{} residual analysis. An $hp$-adaptive procedure and a series of numerical experiments illustrating the performance of the proposed idea will be presented in Section \[sc:adaptivity\]. Finally, we summarize our findings in Section \[sc:concl\].

Adaptive Newton Methods in Banach Spaces {#sc:anewton}
========================================

In the following section we shall briefly revisit an adaptive prediction-type Newton algorithm from [@AmreinWihler:15].

Abstract Framework
------------------

Let $ X, Y $ be two Banach spaces, with norms $\|\cdot\|_X$ and $\|\cdot\|_Y$, respectively. Given an open subset $\Xi\subset X$, and a (possibly nonlinear) operator ${\mathsf{F}}:\,\Xi\to Y$, we consider the [*nonlinear*]{} operator equation $$\label{eq:F=0}
{\mathsf{F}}(u)=0,$$ for some unknown zeros $u\in\Xi$. Supposing that the Fréchet derivative ${\mathsf{F}}'$ of ${\mathsf{F}}$ exists in $\Xi$ (or in a suitable subset), the classical Newton method for solving  starts from an initial guess $u_0\in\Xi$, and generates the iterates $$\label{eq:newton}
u_{n+1}=u_n+\delta_n,\qquad n\ge 0,$$ where the update $\delta_n\in X$ is implicitly given by the [*linear*]{} equation $${\mathsf{F}}'(u_n)\delta_n=-{\mathsf{F}}(u_n),\qquad n\ge 0.$$ Naturally, for this iteration to be well-defined, we need to assume that ${\mathsf{F}}'(u_n)$ is invertible for all $n\ge 0$, and that $\{u_n\}_{n\ge 0}\subset\Xi$.

An Adaptive Prediction Strategy {#sc:newton}
-------------------------------

In order to improve the reliability of the Newton method  in the case that the initial guess $u_0$ is relatively far away from a root $u_{\infty}\in\Xi$ of ${\mathsf{F}}$, ${\mathsf{F}}(u_\infty)=0$, introducing some damping in the Newton method is a well-known remedy. In that case  is rewritten as $$\label{eq:damped}
u_{n+1}=u_n-\Delta t_n{\mathsf{F}}'(u_n)^{-1}{\mathsf{F}}(u_n),\qquad n\ge 0,$$ where $\Delta t_n>0$, $n\ge 0$, is a damping parameter that may be adjusted [*adaptively*]{} in each iteration step.

Provided that ${\mathsf{F}}'(u)$ is invertible on a suitable subset of $\Xi\subset X$, we define the [*Newton transform*]{} by $$u\mapsto{\mathsf{N}_{{\mathsf{F}}}}(u):=-{\mathsf{F}}'(u)^{-1}{\mathsf{F}}(u).$$ Then, rearranging terms in , we notice that $$\frac{u_{n+1}-u_{n}}{\Delta t_n}={\mathsf{N}_{{\mathsf{F}}}}(u_n), \qquad n\ge 0,$$ i.e., can be seen as the discretization of the [*Davydenko-type system*]{}, $$\label{eq:davy}
\begin{split}
\dot{u}(t)&={\mathsf{N}_{{\mathsf{F}}}}(u(t)), \quad t\geq 0,\qquad
 u(0)=u_0,
 \end{split}$$ by the forward Euler scheme, with step size $\Delta t_n>0$.

For $t\in[0,\infty)$, the solution $u(t)$ of , if it exists, defines a trajectory in $X$ that starts at $u_0$, and that will potentially converge to a zero of ${\mathsf{F}}$ as $t\to\infty$. Indeed, this can be seen (formally) from the integral form of , that is, $${\mathsf{F}}(u(t))={\mathsf{F}}(u_0)e^{-t},\qquad t\ge 0,$$ which implies that ${\mathsf{F}}(u(t))\to 0$ as $t\to\infty$.

Now taking the view of dynamical systems, our goal is to compute an upper bound for the value of the step sizes $\Delta t_n>0$ from , $n\ge 0$, so that the discrete forward Euler solution $\{u_n\}_{n\ge 0}$ from  stays reasonably close to the continuous solution of . Specifically, for a prescribed tolerance $\tau>0$, a Taylor expansion analysis (see [@AmreinWihler:15 Section 2] for details) reveals that $$u(t)=u_0+t{\mathsf{N}_{{\mathsf{F}}}}(u_0)+\frac{t^2}{2h}\eta_h+\mathcal{O}(t^3)+\mathcal{O}(t^2h\|{\mathsf{N}_{{\mathsf{F}}}}(u_0)\|_X^2),$$ where, for any sufficiently small $h>0$, we let $\eta_h={\mathsf{N}_{{\mathsf{F}}}}(u_0+h{\mathsf{N}_{{\mathsf{F}}}}(u_0))-{\mathsf{N}_{{\mathsf{F}}}}(u_0)$. Hence, after the first time step of length $\Delta t_0>0$ there holds $$\label{eq:O}
u(\Delta t_0)-u_1=\frac{\Delta t_0^2}{2h}\eta_h+\mathcal{O}(\Delta t_0^3)+\mathcal{O}(\Delta t_0^2h\|{\mathsf{N}_{{\mathsf{F}}}}(u_0)\|_X^2),$$ where $u_1$ is the forward Euler solution from . Therefore, upon setting $$\Delta t_0=\sqrt{2\tau h\|\eta_h\|_X^{-1}},$$ we arrive at $$\|u(\Delta t_0)-u_1\|_X\le\tau+\mathcal{O}(\Delta t_0^3)+\mathcal{O}(\Delta t_0^2h\|{\mathsf{N}_{{\mathsf{F}}}}(u_0)\|_X^2).$$ In order to balance the $\mathcal{O}$-terms in  it is sensible to make the choice $$h=\mathcal{O}(\Delta t_0\|{\mathsf{N}_{{\mathsf{F}}}}(u_0)\|_X^{-2}),$$ i.e., $$\label{eq:h}
h=\gamma \Delta t_0\|{\mathsf{N}_{{\mathsf{F}}}}(u_0)\|_X^{-2},$$ for some parameter $\gamma>0$. This leads to the following *adaptive Newton algorithm*.\

 \[al:zs\] Fix a tolerance $\tau>0$ as well as a parameter $\gamma>0$, and set $n\gets0$.

Start the Newton iteration with an initial guess $ u_{0}\in\Xi$. let $\kappa_{n}=\Delta t_{n-1}$, and $h_n=\gamma \kappa_n\|{\mathsf{N}_{{\mathsf{F}}}}(u_n)\|_X^{-2}$ based on ; define the Newton step size $$\label{eq:knew}
\Delta t_n=\min\left\{\sqrt{2\tau h_n{\left\|{\mathsf{N}_{{\mathsf{F}}}}(u_n+h_n{\mathsf{N}_{{\mathsf{F}}}}(u_n))-{\mathsf{N}_{{\mathsf{F}}}}(u_n)\right\|}^{-1}_{X}},1\right\}.$$ Compute $u_{n+1}$ based on the Newton iteration , and go to ([3:]{}) with $n\leftarrow n+1 $.

The following comments are noteworthy.

(i) We remark that the purpose of the above derivations is to work under minimal structural assumptions on the dynamics caused by the Newton transform. This is particularly important in the context of nonlinear singularly perturbed problems (with ${\epsilon}\ll 1$): Indeed, for such problems (as ${\epsilon}\to 0$) the inverse of the derivative of the associated nonlinear operator will typically become highly irregular (even on a local level), and, for instance, local or global Lipschitz type or boundedness assumptions (involving, e.g., second derivative information) will in general not be available with practically feasible constants. We remark, however, that more sophisticated step size prediction schemes can be derived if the structure of the nonlinearities can be exploited further; see, in particular, the monograph [@5].

(ii) The minimum in  ensures that the step size $\Delta t_n$ is chosen to be 1 whenever possible. Indeed, this is required in order to guarantee quadratic convergence of the Newton iteration close to a root (provided that the root is simple).

(iii) The preset tolerance $ \tau $ in the above adaptive strategy will typically be fixed [*a priori*]{}. Here, for highly nonlinear problems featuring numerous or even infinitely many solutions, it is typically mandatory to select $\tau\ll 1$ small in order to remain within the attractor of the given initial guess. This is particularly important if the starting value is relatively far away from a solution.

    ------------------------------------------------------------------------

Newton-Galerkin $hp$-FEM for Semilinear Elliptic Problems {#sc:FEM}
=========================================================

The purpose of this section is to derive an $hp$-version Newton-Galerkin finite element formulation for the solution of . To this end, we will apply the abstract adaptive Newton framework from the previous Section \[sc:anewton\], and apply an $hp$-discretization methodology for the resulting sequence of (locally) linearized problems.

Newton Linearization
--------------------

In order to apply an adaptive Newton method as introduced in Section \[sc:anewton\] to the nonlinear problem , note that the Fréchet-derivative of $ {\mathsf{F}}_{{\epsilon}} $ from  at $u\in X$ is given by $$\left \langle {\mathsf{F}}_{{\epsilon}}'(u)w,v\right \rangle = \int_{\Omega}\left\{{\epsilon}w'v'-f'(u)wv\right\}{\,{\mathsf{d}}x},\qquad v,w\in X=H^1_0(\Omega),$$ where we write $f'\equiv\partial_uf$. We note that, if $f'(u)\in L^{1}(\Omega)$, then ${\mathsf{F}}_{{\epsilon}}'(u)$ is a well-defined linear and bounded mapping from $X$ to $X'$; see [@AmreinWihler:15 Lemma A.1].

Then, given an initial guess $u_0\in X$ for , the Newton method  is to find $u_{n+1}\in X$ from $u_n\in X$, $n\ge 0$, such that $${\mathsf{F}}_{\epsilon}'(u_n)(u_{n+1}-u_n)=-\Delta t_n{\mathsf{F}}_{\epsilon}(u_n),$$ in $X'$. Equivalently, $$\label{eq:weak}
a_{\epsilon}(u_n; u_{n+1},v)=a_{\epsilon}(u_n; u_{n},v)-\Delta t_n\ell_{\epsilon}(u_n;v)\qquad\forall v\in X,$$ where, for [*fixed*]{} $u\in X$, $$\begin{aligned}
a_{{\epsilon}}(u;w,v)&:=\int_{\Omega}\left\{{\epsilon}w'v'-f'(u)wv\right\}{\,{\mathsf{d}}x},\\
l_{{\epsilon}}(u;v)&:=\int_{\Omega}\left\{{\epsilon}u'v'-f(u)v\right\}{\,{\mathsf{d}}x}\end{aligned}$$ are bilinear and linear forms on $X\times X$ and $X$, respectively.

Incidentally, if there are constants $\underline{\lambda},\overline{\lambda}\ge 0$ with ${\epsilon}C_P^{-2}>\overline{\lambda}$ such that $ -\underline\lambda \leq f'(u) \leq \overline{\lambda}$ holds for all $u\in\mathbb{R} $, where $C_P:=\nicefrac{(b-a)}{\pi}>0$ is the (best) constant in the Poincaré inequality on $\Omega=(a,b)$, i.e.,$$\label{eq:P}
\|w\|_{0}\le C_P\|w'\|_{0}\qquad\forall w \in X,$$ see [@HLP52 Theorem 257], then, it is elementary to show that the formulation  is a linear second-order diffusion-reaction problem, which is coercive on $X$. In consequence, for any given $ u_n \in X $, it has a unique solution $u_{n+1}\in X$; see [@AmreinDiss15] for details.

------------------------------------------------------------------------

Linearized Finite Element Approach
----------------------------------

In order to provide a numerical approximation of , we will discretize the weak formulation  by means of an $hp$-finite element method, which, in combination with the Newton iteration, constitutes a Newton-Galerkin approximation scheme. Furthermore, in view of deriving an $hp$-adaptive algorithm, we will develop an ${\epsilon}$-robust $hp$-version [*a posteriori*]{} residual analysis for the linear finite element formulation.

### $hp$-Spaces

We introduce a partition ${\mathcal{T}}=\{K_j\}_{j=1}^N$ of $N\ge 1$ (open) elements $K_j=(x_{j-1},x_j)$, $j=1,2,\ldots,N$, on $\Omega=(a,b)$, with $a=x_0<x_1<x_2<\ldots<x_{N-1}<x_N=b$. The length of an element $K_j$ is denoted by $h_j=x_j-x_{j-1}$, $j=1,2,\ldots,N$. For each element $K_j \in {\mathcal{T}}$, it will be convenient to introduce the patch $\widetilde K_j = \bigcup \{K'\in{\mathcal{T}}\,|\, 
\overline{K'}\cap \overline{K_j} \ne \emptyset\}$ as the union of $K_j$ and of the elements adjacent to it. In addition, to each element $K_j$ we associate a polynomial degree $p_j\ge 1$, $j=1,2,\ldots,N$. These numbers are stored in a polynomial degree vector ${\bm p}=(p_1,p_2,\ldots,p_N)$. Then, we define an $hp$-finite element space by $${V^{\rm hp}({\mathcal{T}},{\bm p})}=\left\{v\in H^1_0(\Omega):\,
v|_{K_j}\in\mathbb{P}_{p_j}(K_j),\,j=1,2,\ldots,N\right\},$$ where, for $p\ge 1$, we denote by $\mathbb{P}_p$ the space of all polynomials of degree at most $p$. We say that the pair $({\mathcal{T}},{\bm p})$ of a partition ${\mathcal{T}}$ and of a degree vector ${\bm p}$ is $\mu$-shape regular, for some constant $\mu>0$ independent of $j$, if $$\label{eq:gamma-shape-regular}
\mu^{-1} h_{j+1} \leq h_{j} \leq \mu h_{j+1}, 
\qquad  
\mu^{-1} p_{j+1} \leq p_{j} \leq \mu p_{j+1}, 
\qquad j=1,\ldots,N-1,$$ i.e., if both the element sizes and polynomial degrees of neighboring elements are comparable.

### Linear $hp$-FEM

We define the numerical approximation of  in an iterative way as follows: For $n=0$, we let ${V^{\rm hp}({\mathcal{T}},{\bm p})}$ be an initial $hp$-discretization space, and ${u^{\rm hp}}_{0} \in {V^{\rm hp}({\mathcal{T}},{\bm p})}$ a suitable initial guess for the Newton iteration. Furthermore, for $n\ge 0$, and given ${u^{\rm hp}}_n \in {V^{\rm hp}({\mathcal{T}},{\bm p})}$, we find a finite element approximation ${u^{\rm hp}}_{n+1}\in {V^{\rm hp}({\mathcal{T}},{\bm p})}$ of  as the solution of the weak formulation $$\label{eq:fem}
a_{\epsilon}({u^{\rm hp}}_n; {u^{\rm hp}}_{n+1},v)=a_{\epsilon}({u^{\rm hp}}_n; {u^{\rm hp}}_n,v)-\Delta t\ell_{\epsilon}({u^{\rm hp}}_n;v)\qquad\forall v\in {V^{\rm hp}({\mathcal{T}},{\bm p})}.$$ Here, $\Delta t$ takes the role of a (possibly varying) step size parameter in the adaptive Newton scheme as described earlier in Section \[sc:newton\].

[*A Posteriori*]{} Residual Analysis
====================================

Recalling , for given $\Delta t>0$, let us introduce the function $$\label{eq:ut}
u_{n+1}^{(\Delta t,{{\rm hp}})}:={u^{\rm hp}}_{n+1}-(1-\Delta t){u^{\rm hp}}_n,$$ as well as $$\label{eq:ft}
f^{\Delta t}({u^{\rm hp}}_{n+1}):=\Delta tf({u^{\rm hp}}_n)+f'({u^{\rm hp}}_n)(u_{n+1}^{{\rm hp}}-{u^{\rm hp}}_n).$$ Then, rearranging terms, we can rewrite  in the following form: $$\label{start}
 \begin{aligned}
{\epsilon}\int_{\Omega}{ (u_{n+1}^{(\Delta t,{{\rm hp}})})'v'} {\,{\mathsf{d}}x}& = \int_{\Omega}{ f^{\Delta t} (u_{n+1}^{{\rm hp}}) v }{\,{\mathsf{d}}x}\qquad \forall v  \in V_{0}^h.
\end{aligned}$$

The aim of this section is to derive [*a posteriori*]{} residual-based bounds for . In order to measure the discrepancy between the finite element discretization  and the original problem , an obvious quantity to bound is the residual ${\mathsf{F}}_{{\epsilon}}({u^{\rm hp}}_{n+1})$ in $X'$. In order to proceed in this direction, we notice that the adaptively chosen damping parameter $\Delta t$ in the Newton method  will equal 1 sufficiently close to a root of ${\mathsf{F}}_{\epsilon}$. For this reason, as was proposed in [@AmreinWihler:15], we may focus on the ‘shifted’ residual ${\mathsf{F}}_{{\epsilon}}(u_{n+1}^{(\Delta t,{{\rm hp}})})$ in $X'$ instead; indeed, estimating this quantity turns out to be very natural in view of the linearized finite element discretization .

Robust Interpolation Estimates {#sc:int}
------------------------------

For the purpose of deriving an (upper) *a posteriori* residual estimate that is robust with respect to the singular perturbation parameter ${\epsilon}$ (for $0<{\epsilon}\ll 1$) as well as optimally scaled with respect to the local element sizes $h_j$ and polynomial degrees $p_j$, it is crucial to construct an interpolation operator that is *simultaneously* $L^2$- and $H^1$-stable. This will be accomplished in the current section; see Proposition \[pr:intp\] and Corollary \[cr:int\] below.

\[pr:intp\] Let the pair $({\mathcal{T}},{\bm p})$ be $\mu$-shape regular; see . Then, there exists an interpolation operator ${\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}:\,H^1_0(\Omega)\to{V^{\rm hp}({\mathcal{T}},{\bm p})}$ such that, for any $j=1,2,\ldots,N$, and for all $v\in H^1_0(\Omega)$, there holds $$\label{eq:CI}
\begin{split}
{\left\|v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v\right\|}_{L^2(K_j)}&\le C_I\min\left\{{\left\|v\right\|}_{L^2(\widetilde K_j)}, h_jp_j^{-1}{\left\|v'\right\|}_{L^2(\widetilde K_j)}\right\},\\
{\left\|(v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v)'\right\|}_{L^2(K_j)}&\le C_I{\left\|v'\right\|}_{L^2(\widetilde K_j)}.
\end{split}$$ Here, $C_I>0$ is a constant that depends solely on $\mu$; in particular, it is independent of $v$, $\mathcal{T}$, and of $\bm p$.

Let us, without loss of generality, assume that $\Omega=(0,1)$. The result can be shown with the techniques employed in the higher-dimensional case in [@karkulik-melenk15 Cor. [3.7]{}]. In the present, one-dimensional case, a significantly simpler argument can be brought to bear, which allows us to make the paper self-contained. Let $x_{-1} = -h_1$ and $x_{N+1} = 1+h_N$ and $\varphi_i$, $i=0,\ldots,N+1$ be the standard piecewise linear hat functions associated with the nodes $x_i$, $i=-1,\ldots,N+1$. The extra nodes $x_{-1}$ and $x_{N+1}$ define in a natural way the elements $K_{0}$ and $K_{N+1}$. The (open) patches $\omega_i$, $i=0,\ldots,N$, are given by the supports of the functions $\varphi_i$, i.e., $\omega_i = (\operatorname*{supp} \varphi_i)^\circ
 = K_i \cup K_{i+1} \cup \{x_i\}$.

Polynomial approximation (see, e.g., [@MelenkCLEM Proposition [A.2]{}]) gives the existence of an interpolation operator $J_p:L^2(-1,1) \rightarrow {\mathbb P}_p(-1,1)$ that is uniformly (in $p\ge 0$) stable, i.e., $\|J_p v\|_{L^2(-1,1)} \leq C \|v\|_{L^2(-1,1)}$ for all $v \in L^2(-1,1)$ and has the following properties for $v \in H^1(-1,1)$: $$(p+1) \|v - J_p v\|_{L^2(-1,1)} + 
\|(v - J_p v)^\prime\|_{L^2(-1,1)} \leq  C \|v^\prime\|_{L^2(-1,1)}.$$ Furthermore, if $v$ is antisymmetric with respect to the midpoint $x = 0$, then $J_p v$ can be assumed to be antisymmetric as well, i.e., $(J_p v)(0) = 0$ (this follows from studying the antisymmetric part of the original function $J_p v$).

The approximation ${\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v$ is now constructed with the aid of a “partition of unity argument” as described in [@babuska-melenk96 Theorem [2.1]{}]. For $\omega_0$ and $\omega_N$, extend $v$ anti-symmetrically, i.e., $v(x):=-v(-x)$ for $x \in K_{0}$ and $v(x):=-v(1-x)$ for $x \in K_{N+1}$. Then $v$ is defined on each patch $\omega_i$, $i=0,\ldots,N$. For each patch $\omega_i$, let $p^\prime_i:= \min\{p_{i},p_{i+1}\}$ (with the understanding $p_0 = p_1$ and $p_{N+1} = p_N$). The above operator $J_p$ then induces for each patch $\omega_i$ by scaling an operator $J^i:L^2(\omega_i) \rightarrow {\mathcal P}_{p^\prime_i-1}(\omega_i)$ with the following properties: $$\frac{p_i^\prime+1}{h_i}\|v - J^i v\|_{L^2(\omega_i)} + 
\|(v - J^i v)^\prime\|_{L^2(\omega_i)} \leq  C \|v^\prime\|_{L^2(\omega_i)};$$ here, we have exploited the $\mu$-shape regularity of the mesh. We note that $(J^0 v)(0) = 0$ and $(J^N v)(1) = 0$. Also, the operators $J^i$ are uniformly (in the polynomial degree) stable in $L^2(\omega_i)$. The approximation ${\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v$ is now taken to be ${\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v := \sum_{i=0}^N \varphi_i J^i v$. The desired approximation properties follow now from [@babuska-melenk96 Theorem [2.1]{}].

The above proposition implies the following bounds.

\[cr:int\] For $v\in H^1_0(\Omega)$, the interpolant from Proposition \[pr:intp\] satisfies $${\left\|v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v\right\|}^2_{L^2(K_j)}\le
C^2_I\alpha_j{\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon},\widetilde K_j}^2 ,
\qquad j=1,2,\ldots,N,$$ and $$\big|(v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v)(x_j)\big|^2\le
C_I^2\gamma_j\left({\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon},\widetilde K_{j}}^2+{\left|\!\left|\!\left|v\right|\!\right|\!\right|}^2_{{\epsilon},\widetilde K_{j+1}}\right),
\qquad j=1,2,\ldots,N-1.$$ Here, we let $$\begin{aligned}
\label{eq:alpha}
\alpha_j&=
\min\left\{{\epsilon}^{-1}h_j^2p_j^{-2},1\right\},\qquad j=1,\ldots, N,\end{aligned}$$ and, with $$\begin{aligned}
\beta_j&=\alpha_jh_j^{-1}+2\sqrt{{\epsilon}^{-1}\alpha_j},\label{eq:beta}\end{aligned}$$ we define $$\label{eq:gamma}
\gamma_j=\frac{\beta_j\beta_{j+1}}{\beta_j+\beta_{j+1}},$$ for $1\le j\le N-1$, and $\gamma_0=\gamma_N=0$. Moreover, $C_I$ is the constant from .

We proceed along the lines of [@Ve1998]. Using the bounds from Proposition \[pr:intp\], for each element $K_j\in{\mathcal{T}}$, we have that $$\begin{aligned}
{\left\|v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v\right\|}^2_{L^2(K_j)}&\le 
C^2_I\frac{h_j^2}{{\epsilon}p_j^2}{\epsilon}{\left\|v'\right\|}^2_{L^2(\widetilde K_j)}.\end{aligned}$$ Furthermore, $$\begin{aligned}
{\left\|v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v\right\|}^2_{L^2(K_j)}
&\le C^2_I{\left\|v\right\|}_{L^2(\widetilde K_j)}^2.\end{aligned}$$ Combining these two estimates, yields the first bound.

In order to prove the second estimate, we apply, for $1\le j\le N-1$, a multiplicative trace inequality; see [@MelenkWihler:15 Lemma A.1]): $$\begin{aligned}
\left|(v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v)(x_j)\right|^2
&\le h_j^{-1}{\left\|v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v\right\|}^2_{L^2(K_j)}\\
&\quad+2{\left\|v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v\right\|}_{L^2(K_j)}{\left\|(v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v)'\right\|}_{L^2(K_j)}.\end{aligned}$$ Then, invoking the above bounds as well as the estimates from Proposition \[pr:intp\], we get $$\begin{aligned}
\big|(v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v)(x_j)\big|^2
&\le C_I^2\left(\alpha_jh_j^{-1}{\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon},\widetilde K_j}^2
+2\sqrt{\alpha_j}{\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon},\widetilde K_j}{\left\|v'\right\|}_{L^2(\widetilde K_j)}\right)\\
&\le C_I^2\left(\alpha_jh_j^{-1}{\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon},\widetilde K_j}^2
+2\sqrt{{\epsilon}^{-1}\alpha_j}{\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon},\widetilde K_j}^2\right)\\
&\le C_I^2\beta_j{\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon},\widetilde K_j}^2,\end{aligned}$$ with $\beta_j$ from . Since $x_j$ is also a boundary point of $K_{j+1}$, we similarly obtain that $$\big|(v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v)(x_j)\big|^2
\le C_I^2\beta_{j+1}{\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon},\widetilde K_{j+1}}^2.$$ Therefore, $$\begin{aligned}
\big|(v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v)(x_j)\big|^2
&=\frac{\beta_{j+1}}{\beta_{j}+\beta_{j+1}}\big|(v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v)(x_j)\big|^2\\
&\quad +\frac{\beta_{j}}{\beta_{j}+\beta_{j+1}}\big|(v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v)(x_j)\big|^2\\ &\le
C_I^2\gamma_j\left({\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon},\widetilde K_{j}}^2+{\left|\!\left|\!\left|v\right|\!\right|\!\right|}^2_{{\epsilon},\widetilde K_{j+1}}\right),\end{aligned}$$ with $\gamma_j$ from . Thus, we have shown the second estimate.

Robust *A Posteriori* Residual Estimate
---------------------------------------

Based on the previous interpolation analysis in Section \[sc:int\], we are now in a position to prove the following main result.

\[thm:main\] \[thm:1\] For $u_{n+1}^{(\Delta t,{{\rm hp}})}$ from  there holds the upper *a posteriori* bound $$\label{eq:upperbound}
{\left|\!\left|\!\left|{\mathsf{F}}_{\epsilon}(u_{n+1}^{(\Delta t,{{\rm hp}})})\right|\!\right|\!\right|}_{X',{\epsilon}}^2
\preccurlyeq \delta_{n,\Omega}^2+\sum_{j=1}^N{\eta_{n,j}^{2}},$$ with $$\label{Femerror}
\begin{split}
\eta_{n,j}^2:&= \alpha_{j} {\left\|f^{\Delta t}({u^{\rm hp}}_{n+1})+{\epsilon}(u_{n+1}^{(\Delta t,{{\rm hp}})})''\right\|}_{0,K_j}^2\\
&\quad+\frac{1}{2}{\epsilon}^2\gamma_j\left|{\left\llbracket(u_{n+1}^{(\Delta t,{{\rm hp}})})'\right\rrbracket}(x_j)\right|^2
+\frac{1}{2}{\epsilon}^2\gamma_{j-1}\left|{\left\llbracket(u_{n+1}^{(\Delta t,{{\rm hp}})})'\right\rrbracket}(x_{j-1})\right|^2,
\end{split}$$ for $1\le j\le N$, and $$\label{Newtonerror}
\delta_{n,\Omega}:={\left\|f^{\Delta t}(u_{n+1}^{{{\rm hp}}})-f(u_{n+1}^{(\Delta t,{{\rm hp}})})\right\|}_{0,\Omega}.$$ Here, $\alpha_j$ and $\gamma_j$ are defined in  and , respectively, and we set $\gamma_0=\gamma_N=0$.

We emphasize that the constants $\alpha_j$ and ${\epsilon}^2\alpha_j\gamma_j$ appearing in the error indicators $\eta_{K_j}$ from  remain bounded as $h_j,{\epsilon}\to 0$ (and $p_j\to\infty$). In addition, we note that $\nicefrac{1}{2} \min\{\beta_j,\beta_{j+1}\} \leq \gamma_j \leq \min\{\beta_j,\beta_{j+1}\}$, and that $2\sqrt{{\epsilon}^{-1} \alpha_j } \leq \beta_j \leq 3 \sqrt{{\epsilon}^{-1} \alpha_j}$.

------------------------------------------------------------------------

For *linear problems*, the (adaptive) Newton linearization is redundant. In particular, upon setting $\Delta t=1$, we infer that $\delta_{n,\Omega}=0$ in , and the residual is bounded by the local error indicators $\eta_{n,j}$ only. These quantities, in turn, can be exploited for the purpose of designing an $hp$-adaptive mesh refinement strategy, as outlined, for example, in Section \[sc:hpadapt\] below. We emphasize that this approach has been studied earlier in the report [@MelenkWihler:15], where a number of numerical tests in the context of $hp$-adaptive FEM for singularly perturbed linear problems have been performed; in particular, these experiments have underlined the robustness of the bound  as ${\epsilon}\to 0$ in the linear case. In Section \[sc:adaptivity\], we will extend this idea by proposing a suitable interplay between $hp$-mesh refinements and the (adaptive) Newton linearization from Algorithm \[al:zs\].

------------------------------------------------------------------------

Given $u_{n+1}^{(\Delta t,{{\rm hp}})}$ and $f^{\Delta t}({u^{\rm hp}}_{n+1})$ from  and , respectively, and any $v\in H^1_0(\Omega)$, there holds the identity $$\label{eq:Fres}
\left \langle {\mathsf{F}}_{{\epsilon}}(u_{n+1}^{(\Delta t,{{\rm hp}})}),v \right \rangle
=\sum_{j=1}^{N-1}a_j+\sum_{j=1}^N(b_j+c_j),$$ where $$\begin{aligned}
{2}
a_j&:= {\epsilon}{\left\llbracket(u_{n+1}^{(\Delta t,{{\rm hp}})})'\right\rrbracket}(x_j)(v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v)(x_j),&\qquad&1\le j\le N-1,\nonumber\\
b_j&:=\int_{K_j}{\left\{f^{\Delta t}({u^{\rm hp}}_{n+1})+{\epsilon}(u_{n+1}^{(\Delta t,{{\rm hp}})})''\right\}({\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v-v)}{\,{\mathsf{d}}x},&&1\le j\le N,\label{bounds}\\
c_j&:=\int_{K_j}{\left\{f^{\Delta t}({u^{\rm hp}}_{n+1})-f(u_{n+1}^{(\Delta t,{{\rm hp}})})\right\}v}{\,{\mathsf{d}}x},&&1\le j\le N;\nonumber\end{aligned}$$ see [@AmreinWihler:15 Proposition 4.3].

We estimate the terms $a_j$, $b_j$, and $c_j$ individually. First let $1\le j\le N-1$. Then, $a_j$ from  can be estimated using Corollary \[cr:int\] as follows: $$\begin{aligned}
|a_j|&\le {\epsilon}\left|{\left\llbracket(u_{n+1}^{(\Delta t,{{\rm hp}})})'\right\rrbracket}(x_j)\right||(v-{\pi_{{V^{\rm hp}({\mathcal{T}},{\bm p})}}}v)(x_j)|\\
&\preccurlyeq\gamma_j^{\nicefrac12}{\epsilon}\left|{\left\llbracket(u_{n+1}^{(\Delta t,{{\rm hp}})})'\right\rrbracket}(x_j)\right|\left({\left|\!\left|\!\left|v\right|\!\right|\!\right|}^2_{{\epsilon},\widetilde K_j}+{\left|\!\left|\!\left|v\right|\!\right|\!\right|}^2_{{\epsilon},\widetilde{K}_{j+1}}\right)^{\nicefrac12}.\end{aligned}$$ Applying the Cauchy-Schwarz inequality leads to $$\begin{aligned}
\left|\sum_{j=1}^{N-1}a_j\right|
&\preccurlyeq\left(\sum_{j=1}^{N-1}{{\epsilon}^2\gamma_j\left|{\left\llbracket(u_{n+1}^{(\Delta t,{{\rm hp}})})'\right\rrbracket}(x_j)\right|^2}\right)^{\nicefrac12}
\left(\sum_{j=1}^{N-1}\left({\left|\!\left|\!\left|v\right|\!\right|\!\right|}^2_{{\epsilon},\widetilde K_j}+{\left|\!\left|\!\left|v\right|\!\right|\!\right|}^2_{{\epsilon},\widetilde{K}_{j+1}}\right)\right)^{\nicefrac12}\\
&\preccurlyeq\left(\sum_{j=1}^{N-1}{{\epsilon}^2\gamma_j\left|{\left\llbracket(u_{n+1}^{(\Delta t,{{\rm hp}})})'\right\rrbracket}(x_j)\right|^2}\right)^{\nicefrac12}{\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon}}.\end{aligned}$$ Furthermore, again using Corollary \[cr:int\], we see that $$\begin{aligned}
\left|\sum_{j=1}^Nb_j\right|
&\preccurlyeq\sum_{j=0}^N\alpha_j^{\nicefrac12}{\left\|f^{\Delta t}({u^{\rm hp}}_{n+1})+{\epsilon}(u_{n+1}^{(\Delta t,{{\rm hp}})})''\right\|}_{0,K_j}{\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon},\widetilde K_j}\\
&\preccurlyeq\left(\sum_{j=1}^N\alpha_j{\left\|f^{\Delta t}({u^{\rm hp}}_{n+1})+{\epsilon}(u_{n+1}^{(\Delta t,{{\rm hp}})})''\right\|}_{0,K_j}^2\right)^{\nicefrac12}\left(\sum_{j=1}^N{\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon},\widetilde{K}_j}^2\right)^{\nicefrac12}\\
&\preccurlyeq\left(\sum_{j=1}^N\alpha_j{\left\|f^{\Delta t}({u^{\rm hp}}_{n+1})+{\epsilon}(u_{n+1}^{(\Delta t,{{\rm hp}})})''\right\|}_{0,K_j}^2\right)^{\nicefrac12}{\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon}}.\end{aligned}$$ Similarly, there holds $$\begin{aligned}
\left|\sum_{j=1}^Nc_j\right|&
\preccurlyeq\sum_{j=1}^N {\left\|f^{\Delta t}({u^{\rm hp}}_{n+1})-f(u_{n+1}^{(\Delta t,{{\rm hp}})})\right\|}_{0,K_j}{\left\|v\right\|}_{0,K_j}\\
&\preccurlyeq\left(\sum_{j=1}^N{\left\|f^{\Delta t}({u^{\rm hp}}_{n+1})-f(u_{n+1}^{(\Delta t,{{\rm hp}})})\right\|}_{0,K_j}^2\right)^{\nicefrac12}{\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon}}.\end{aligned}$$ Now, applying the Cauchy-Schwarz inequality to  we see that $$\begin{aligned}
{\left|\left \langle {\mathsf{F}}_{{\epsilon}}(u_{n+1}^{(\Delta t,{{\rm hp}})}),v \right \rangle\right|}
&\preccurlyeq\left(\delta_{n,\Omega}^2+\sum_{j=1}^N{\eta_{n,j}^{2}}\right)^{\nicefrac{1}{2}}{\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon}}.
\end{aligned}$$ Dividing by ${\left|\!\left|\!\left|v\right|\!\right|\!\right|}_{{\epsilon}}$, and taking the supremum for all $v\in H^1_0(\Omega)$, completes the proof.

Suppose that there exist constants $\lambda>-C_P^{-2}{\epsilon}$ and $L>0$, with $C_P$ being the Poincaré constant from , such that $$\begin{aligned}
(f(x)-f(y))(x-y)\le-\lambda(x-y)^2,\qquad
|f(x)-f(y)|\le L|x-y|,\end{aligned}$$ for all $x,y\in\mathbb{R}$. Then, the residual norm $ {\left|\!\left|\!\left|{\mathsf{F}}_{\epsilon}(u_{n+1}^{(\Delta t,h)})\right|\!\right|\!\right|}_{X',{\epsilon}}$ and the energy norm ${\left|\!\left|\!\left|u-u_{n+1}^{(\Delta t,h)}\right|\!\right|\!\right|}_{{\epsilon}}$ of the error are equivalent; see [@AmreinWihler:15 Theorem 4.5] for details.

------------------------------------------------------------------------

Following along the lines of [@AmreinWihler:15 §4.4.2] and [@MelenkWihler:15 Appendix], it is possible to prove *${\epsilon}$-robust local lower residual bounds* in terms of the error indicators $\eta_{K_j}$ and the data oscillation terms. This approach, however, results in local efficiency bounds that will be slightly suboptimal with respect to the local polynomial degrees due to the need of applying $p$-dependent norm equivalences (involving cut-off functions).

------------------------------------------------------------------------

An $hp$-Adaptive Newton-Galerkin Procedure {#sc:adaptivity}
==========================================

In this section, we will discuss how the *a posteriori* bound from Theorem \[thm:main\] can be exploited for the purpose of designing an $hp$-adaptive Newton-Galerkin algorithm for the numerical solution of . We will begin by revisiting an $hp$-adaptive testing strategy from [@FaWiWi14] (see also [@Wi11; @Wi11_2]).

$hp$-Adaptive Refinement Procedure {#sc:hpadapt}
----------------------------------

In order to $hp$-refine the $hp$-finite element space ${V^{\rm hp}({\mathcal{T}},{\bm p})}$, we shall apply an $hp$-adaptive algorithm which is based on the following two ingredients.

### (a) Element marking: {#a-element-marking .unnumbered}

The elementwise error indicators $\eta_{n,j}$ from Theorem \[thm:main\] are employed in order to mark elements for refinement. More precisely, we fix a parameter $\vartheta\in(0,1)$, and select elements to be refined according to the *Dörfler marking* criterion: $$\label{eq:dorf}\tag{D}
\vartheta\sum_{j=1}^N\eta_{n,j}^2\le \sum_{j'=1}^M\eta_{n,j'}^2.$$ Here, the indices $j'$ are chosen such that the error indicators $\eta_{j'}$ from  are sorted in descending order, and $M$ is minimal.

### (b) $hp$-refinement criterion: {#b-hp-refinement-criterion .unnumbered}

The decision of whether a marked element in step (a) is refined with respect to $h$ (element bisection) or $p$ (increasing the local polynomial order by 1) is based on a smoothness testing approach. Specifically, if the (numerical) solution is considered smooth on a marked element $K_j$, then the polynomial degree is increased by 1 on that particular element (no element bisection), otherwise the element is bisected (retaining the current polynomial degree $p_j$ on both subelements). In order to evaluate the smoothness of the solution ${u^{(\Delta t,\rm hp)}_{n+1}}$ on a marked element $K_j$, we employ an elementwise smoothness indicator as introduced in [@FaWiWi14 Eq. (3)], $$\label{eq:hp}\tag{F}
\mathcal{F}^{p_j}_j\left[{u^{(\Delta t,\rm hp)}_{n+1}}\right]:=
\frac{{\left\|\frac{{\mathsf{d}}^{p_j-1}}{{\mathsf{d}}x^{p_j-1}}{u^{(\Delta t,\rm hp)}_{n+1}}\right\|}_{L^\infty(K_j)}}{h_j^{-\nicefrac12}{\left\|\frac{{\mathsf{d}}^{p_j-1} {u^{(\Delta t,\rm hp)}_{n+1}}}{{\mathsf{d}}x^{p_j-1}}\right\|}_{L^2(K_j)}+\frac{1}{\sqrt2}h_j^{\nicefrac12}{\left\|\frac{{\mathsf{d}}^{p_j} {u^{(\Delta t,\rm hp)}_{n+1}}}{{\mathsf{d}}x^{p_j}}\right\|}_{L^2(K_j)}},$$ if $\frac{{\mathsf{d}}^{p_j-1}}{{\mathsf{d}}x^{p-1}}{u^{(\Delta t,\rm hp)}_{n+1}}|_{K_j}\not\equiv 0$, and $\mathcal{F}^{p_j}_j[{u^{(\Delta t,\rm hp)}_{n+1}}]=1$ otherwise. Here, the basic idea is to consider the continuous Sobolev embedding $H^1(K_j)\hookrightarrow L^\infty(K_j)$, which implies that $$\sup_{v\in H^1(K_j)}\frac{{\left\|v\right\|}_{L^\infty(K_j)}}{h_j^{-\nicefrac12}{\left\|v\right\|}_{L^2(K_j)}+\frac{1}{\sqrt2}h_j^{\nicefrac12}{\left\|v'\right\|}_{L^2(K_j)}}\le 1;$$ see [@FaWiWi14 Proposition 1]. In particular, it follows that $\mathcal{F}_j^{p_j}[{u^{(\Delta t,\rm hp)}_{n+1}}]\le 1$. For ease of evaluation, note that, by taking the derivative of order $p_j-1$ in the definition , the smoothness indicator $\mathcal{F}^{p_j}_j[{u^{(\Delta t,\rm hp)}_{n+1}}]$ is evaluated for linear functions only; in this case, it can be shown that $$\frac12\approx\frac{\sqrt{3}}{\sqrt{6}+1}\le\mathcal{F}_j^{p_j}\left[{u^{(\Delta t,\rm hp)}_{n+1}}\right]\le 1;$$ cf. [@FaWiWi14 Section 2.2]. For a given smoothness testing parameter $\zeta\in(\nicefrac{\sqrt{3}}{(\sqrt{6}+1)},1)$, the numerical solution ${u^{(\Delta t,\rm hp)}_{n+1}}$ is classified smooth on $K_j$ if $$\label{eq:st}
\mathcal{F}^{p_j}_j\left[{u^{(\Delta t,\rm hp)}_{n+1}}\right]\ge\zeta,$$ and otherwise nonsmooth. Incidentally, representing the local solution ${u^{(\Delta t,\rm hp)}_{n+1}}|_{K_j}$ in terms of (local) Legendre polynomials (or more general Jacobi polynomials), any derivatives of ${u^{(\Delta t,\rm hp)}_{n+1}}$ can be evaluated exactly by means of appropriate recurrence relations. For further information, we refer to the paper [@FaWiWi14] (see also [@Wi11_2; @Wi11]), where the smoothness testing strategy described above has been introduced and discussed in detail.

Fully Adaptive Newton-Galerkin Method
-------------------------------------

We will now propose a procedure that provides an interplay of the adaptive Newton methods presented in Section \[sc:anewton\] and automatic $hp$-finite element mesh refinements based on the [*a posteriori*]{} error estimate from Theorem \[thm:1\] (as outlined in the previous Section \[sc:hpadapt\]). To this end, we make the assumption that the Newton sequence $\left\{u_{n+1}^{(\Delta t_n,{{\rm hp}})}\right\}_{n\ge 0} $ given by  and , with a (possibly varying) step size $\Delta t_n$, is well-defined as long as the iterations are being performed.

\[al:full\] Given a (coarse) starting mesh $\mathcal{T}$ in $\Omega$, with an associated (low-order) polynomial degree distribution $\bm p$, and an initial guess $ u_{0}^{{{\rm hp}}} \in  {V^{\rm hp}({\mathcal{T}},{\bm p})}$. Set $n\gets 0$.

Determine the Newton step size parameter $\Delta t_n$ based on ${u^{\rm hp}}_n$ by the adaptive procedure from Algorithm \[al:zs\]. Compute the FEM solution $u_{n+1}^{{{\rm hp}}}$ from  with step size $\Delta t_n$ on the space ${V^{\rm hp}({\mathcal{T}},{\bm p})}$. Furthermore, obtain $ u_{n+1}^{(\Delta t_n,{{\rm hp}})} $ in , and evaluate the corresponding error indicators $ \{\eta_{n,j}\}_j $, and $\delta_{n,\Omega} $ from  and , respectively.

Numerical Experiments
---------------------

Let us provide some comments on how Algorithm \[al:full\] is implemented in the ensuing examples.

-   The Newton transform ${\mathsf{N}_{{\mathsf{F}}}}({u^{\rm hp}}_n)$ required for the computation of the step size parameter $\Delta t_n$ in step ([1:]{}) is approximated using the $hp$-FEM on the current mesh.

-   If the criterion  in step ([3:]{}) is satisfied, the mesh is refined, and then step ([2:]{}) is repeated based on the previously computed solution $u_{n+1}^{{{\rm hp}}}$ as interpolated on the refined mesh.

-   The linear systems resulting from the finite element discretization  are solved by means of a direct solver; in this way, this approach differs from inexact Newton methods as discussed, e.g., in [@ErVo13]; see also [@5].

-   The parameters are chosen to be $\gamma=0.5$ and $\tau=0.1$ in Algorithm \[al:zs\], $\vartheta=0.5$ in , $\zeta=0.6$ in , and $\theta=0.5$ in . Unless stated differently, the procedure is initiated with a uniform initial mesh ${\mathcal{T}}$ consisting of 10 elements, and a polynomial degree distribution $\bm p=(1,\ldots 1)$.

\[ex:slin0\] We begin by looking at the (non-singularly perturbed) *Bratu* equation given by $$\label{eq:Bratu}
-u''=\exp(u+1)\quad\text{on }\Omega=(0,1),\qquad u(0)=u(1)=0.$$ This problem features two positive and concave solutions, which we aim to find by means of the fully adaptive Newton-Galerkin procedure introduced above (cf. also [@AmreinWihler:14 Example 3.4]). To this end, we let the initial guesses to be ${u^{\rm hp}}_0\equiv0$ and ${u^{\rm hp}}_0(x)=\pi_1(10x(1-x))$, where $\pi_1$ denotes the piecewise linear interpolation operator on the initial mesh. For these choices, Algorithm \[al:full\] generates the two solutions in bold and dashed lines, respectively, in Figure \[fig:slin0\] (left). Since both solutions are smooth, the $hp$-adaptive testing strategy  dictates pure $p$-refinements in all steps as expected, and, consequently, the residuals decay at an exponential rate with respect to the number of degrees of freedom; see Figure \[fig:slin0\] (right).

------------------------------------------------------------------------

![Example \[ex:slin0\]: Two positive solutions of  (left), and the corresponding performance data for the residuals (right).[]{data-label="fig:slin0"}](figures/Bratusolutions "fig:"){width="0.45\linewidth"} ![Example \[ex:slin0\]: Two positive solutions of  (left), and the corresponding performance data for the residuals (right).[]{data-label="fig:slin0"}](figures/Braturesiduals "fig:"){width="0.45\linewidth" height="0.375\linewidth"}

\[ex:slin1\] The focus of the next sequence of experiments is on the robustness of the *a posteriori* residual bound  with respect to the singular perturbation parameter ${\epsilon}$ as ${\epsilon}\to 0$. To this end, let us consider the *Ginzburg-Landau*-type equation, $$-{\epsilon}u''=u-u^3\ \text{on} \ (0,1),\qquad u(0)=u(1)=0.$$ Starting from the initial guess depicted in Figure \[fig:slin1\_1\] (left), we test the fully adaptive Newton-Galerkin Algorithm \[al:full\] for different choices of ${\epsilon}=\{10^{-i}\}_{i=2}^5$. As ${\epsilon}\to 0$ the resulting solutions feature ever stronger boundary layers, as well as an interior shock close to the center of the domain; see Figure \[fig:slin1\_1\] (right). The performance data in Figure \[fig:slin1\_2\] (left) shows that the residuals decay, firstly, fairly robust in ${\epsilon}$, and, secondly, exponentially fast with respect to the (square root) of the number of degrees of freedom (cf., e.g., [@Schwab98]). The final mesh for ${\epsilon}=10^{-5}$ is displayed in Figure \[fig:slin1\_2\] (right); whilst the singular effects have been properly resolved by suitable $hp$-refined elements, we notice that the initial finite element space remains nearly unrefined in areas where the solution is approximately constant $\pm 1$.

------------------------------------------------------------------------

![Example \[ex:slin1\] for ${\epsilon}=10^{-5}$: Piecewise linear initial guess (left), and numerical solution with strong layers at the boundary and close to the center of the domain (right).[]{data-label="fig:slin1_1"}](figures/GLinitial "fig:"){width="0.45\linewidth"} ![Example \[ex:slin1\] for ${\epsilon}=10^{-5}$: Piecewise linear initial guess (left), and numerical solution with strong layers at the boundary and close to the center of the domain (right).[]{data-label="fig:slin1_1"}](figures/GLsolution5 "fig:"){width="0.45\linewidth"}

![Example \[ex:slin1\]: Estimated residuals for different choices of ${\epsilon}$ (left), and adaptively generated final $hp$-mesh for ${\epsilon}=10^{-5}$ (right).[]{data-label="fig:slin1_2"}](figures/GLresidual "fig:"){width="0.45\linewidth"} ![Example \[ex:slin1\]: Estimated residuals for different choices of ${\epsilon}$ (left), and adaptively generated final $hp$-mesh for ${\epsilon}=10^{-5}$ (right).[]{data-label="fig:slin1_2"}](figures/GLhpmesh5 "fig:"){width="0.45\linewidth" height="0.355\linewidth"}

\[ex:slin2\] Finally, we turn to the *Fisher* equation, $$\begin{aligned}
-{\epsilon}u'' =u-u^2 \ \text{on} \ (0,1),\qquad
u_{{\epsilon}}(0)=\alpha,\quad \ u_{{\epsilon}}(1)=\beta. 
\end{aligned}$$ Solutions of this problem possess boundary layers, and, furthermore, for $ \alpha>-\nicefrac12 $ and $\beta < 1 $, the solutions feature an increasing number of spikes (which are bounded between 0 and 1) as ${\epsilon}\to 0$; see, e.g., [@Verhulst] for a more detailed discussion. The purpose of this example is to show that our fully adaptive Newton-Galerkin procedure is able to transport appropriate initial guesses sufficiently close along the corresponding trajectory of the Davydenko system , and, thereby, to obtain numerical solutions that retain the topological structure of the starting function. To proceed in this direction, we look at the two initial guesses displayed in Figures \[fig:slin2\_1\] and \[fig:slin2\_2\] (left), as well as at the corresponding numerical solutions (right) resulting from Algorithm \[al:full\]. The initial mesh is based on 20 uniform elements in these experiments. We clearly see that the proposed method has generated two different solutions, which, essentially, show the same *qualitative* properties as the respective initial guesses. We emphasize that, for ${\epsilon}\ll 1$, this will typically not work if a standard Newton approach is employed. Moreover, for the solution given in Figure \[fig:slin2\_2\], the $hp$-refined mesh is plotted in Figure \[fig:slin2\_3\] (left), and shows that the boundary layers and interior spikes have been carefully resolved by the $hp$-adaptive scheme. In addition, from Figure \[fig:slin2\_3\] (right), exponential convergence can be observed for both cases.

------------------------------------------------------------------------

![Example \[ex:slin2\] for ${\epsilon}=0.00025$: Piecewise linear initial guess (left), and numerical solution (right).[]{data-label="fig:slin2_1"}](figures/Finitial1 "fig:"){width="0.45\linewidth"} ![Example \[ex:slin2\] for ${\epsilon}=0.00025$: Piecewise linear initial guess (left), and numerical solution (right).[]{data-label="fig:slin2_1"}](figures/Fsolution1 "fig:"){width="0.45\linewidth"}

![Example \[ex:slin2\] for ${\epsilon}=0.00025$: Piecewise linear initial guess (left), and numerical solution (right).[]{data-label="fig:slin2_2"}](figures/Finitial2 "fig:"){width="0.45\linewidth"} ![Example \[ex:slin2\] for ${\epsilon}=0.00025$: Piecewise linear initial guess (left), and numerical solution (right).[]{data-label="fig:slin2_2"}](figures/Fsolution2 "fig:"){width="0.45\linewidth"}

![Example \[ex:slin2\] for ${\epsilon}=0.00025$: Final $hp$-refined mesh (left), and exponential decay of residuals for the two solutions from Figure \[fig:slin2\_1\] and Figure \[fig:slin2\_2\] (right).[]{data-label="fig:slin2_3"}](figures/Fhpmesh2 "fig:"){width="0.45\linewidth"} ![Example \[ex:slin2\] for ${\epsilon}=0.00025$: Final $hp$-refined mesh (left), and exponential decay of residuals for the two solutions from Figure \[fig:slin2\_1\] and Figure \[fig:slin2\_2\] (right).[]{data-label="fig:slin2_3"}](figures/Fresidual "fig:"){width="0.45\linewidth" height="0.365\linewidth"}

Conclusions {#sc:concl}
===========

The aim of this paper was to introduce a reliable and computationally feasible procedure for the numerical solution of general, semilinear elliptic boundary value problems with possible singular perturbations. The key idea is to combine an adaptive Newton method with an automatic mesh refinement finite element procedure. Here, the (local) Newton damping parameter is selected based on interpreting the scheme within the context of step size control for dynamical systems. Furthermore, the sequence of linear problems resulting from the Newton discretization is treated by means of a robust (with respect to the singular perturbations) $hp$-version [*a posteriori*]{} residual analysis, and by a corresponding adaptive $hp$-refinement scheme. The numerical experiments clearly illustrate the ability of our approach to reliably find solutions reasonably close to the initial guesses (in the sense that solutions with similar qualitative structure like the initial guesses can be found), and to robustly resolve the singular perturbations at an exponential rate.
