1
%\chapter{The $p$-laplacian problem}
1
%\chapter{The $p$-Laplacian problem}
3
3
% ----------------------------------
4
4
\section{Problem statement}
5
5
% ----------------------------------
6
6
\label{sec-p-laplacian}%
8
Let us consider the classical $p$-laplacian problem
8
Let us consider the classical $p$-Laplacian problem
9
9
with homogeneous Dirichlet boundary conditions
10
10
in a domain bounded $\Omega \subset R^d$,
19
19
\bcindex{Dirichlet}%
21
21
When $p=2$, this problem reduces to the linear Poisson problem
22
with homogeneous Dichlet boundary conditions.
23
Otherwise, for any $p>1$, the nonlinear problem is equivalent to the following minimisation problem:
22
with homogeneous Dirichlet boundary conditions.
23
Otherwise, for any $p>1$, the nonlinear problem is equivalent to the following minimization problem:
25
25
{\it (MP): find $u\in W^{1,p}_0(\Omega)$ such that:}
68
68
\section{The fixed-point algorithm}
69
69
% ----------------------------------
70
70
\subsection{Principe of the algorithm}
71
\cindex{algorithm!fixed-point}%
71
\cindex{method!fixed-point}%
72
72
This nonlinear problem is then reduced to a sequence
73
73
of linear subproblems by using the fixed-point algorithm.
74
74
The sequence $\left(u^{(n)}\right)_{n\geq 0}$ is defined
93
93
Let us introduce a mesh ${\cal T}_h$ of $\Omega$
94
94
and the finite dimensional space $X_h$ of continuous
95
picewise polynomial functions
95
piecewise polynomial functions
96
96
and $V_h$, the subspace of $X_h$
97
97
containing elements that vanishes on the boundary of $\Omega$:
110
110
a\left(u^{(n)}_h ;u^{(n+1)}_h,v_h\right) = m(f,v_h),
111
111
\ \forall v_h \in V_h
113
By developping $u_h$ on a basis of $V_h$, this
113
By developing $u_h$ on a basis of $V_h$, this
114
114
problem reduces to a linear system.
115
115
\cindex{form!weighted}%
116
116
The implementation with \Rheolef, involving weighted
135
135
\subsection{Comments}
136
136
\pbindex{Poisson}%
137
The fixed-point algorithm is amorced with $u^{(0)}$ as
137
The fixed-point algorithm is initiated with $u^{(0)}$ as
138
138
the solution of the linear problem associated to $p=2$,
139
139
i.e. the standard Poisson problem with Dirichlet boundary
170
170
\includegraphics[height=6.0cm]{p-laplacian-square-p=1,2-cut.pdf}
173
\caption{The $p$-laplacian for $d=2$:
173
\caption{The $p$-Laplacian for $d=2$:
174
174
(a) elevation view for $p=1.2$;
175
(b) cut along the first bissectice $x_0-x_1=0$.}
175
(b) cut along the first bisector $x_0-x_1=0$.}
176
176
\label{fig-p-laplacian}
178
178
Run the field visualization:
187
187
\cindex{visualization!elevation view}%
188
188
The first command shows an elevation view of the
189
189
solution (see~\ref{fig-p-laplacian}.a) while the second one
190
shows a cut along the first bissectrice $x_0=x_1$.
190
shows a cut along the first bisector $x_0=x_1$.
191
191
(see~\ref{fig-p-laplacian}.b).
193
193
\subsection{Convergence properties of the fixed-point algorithm}
207
207
\includegraphics{p-laplacian-rate-log.pdf}
210
\caption{The fixed-point algorithm on the $p$-laplacian for $d=2$:
210
\caption{The fixed-point algorithm on the $p$-Laplacian for $d=2$:
211
211
(a) convergence when $p<2$;
213
213
(c) convergence rate versus $p$;
268
268
where xdof$(V_h)$ denotes the set of nodes associated to the $V_h$
269
269
degrees of freedom. Since elements of $V_h$ vanishes on the
270
270
boundary, the xdof$(V_h)$ contains all nodes associated to the
271
degrees of freedoms of $X_h$ except nodes located on the boudnary.
271
degrees of freedoms of $X_h$ except nodes located on the boundary.
272
272
Fig~\ref{fig-p-laplacian-rate}.a and~\ref{fig-p-laplacian-rate}.b
273
273
show that the residual term decreases exponentially versus $n$,
274
274
since the slope of the plot in semi-log scale tends to be strait.
299
299
\section{The Newton algorithm}
300
300
% ----------------------------------
301
301
\subsection{Principe of the algorithm}
302
\cindex{algorithm!Newton}%
302
\cindex{method!Newton}%
303
303
An alternative to the fixed-point algorithm is
304
304
to solve the nonlinear problem $(P)$
305
305
by using the Newton algorithm.
498
498
in the file \file{newton.h}.
499
499
The main program is~\file{p\_laplacian\_newton.cc}, that
500
500
uses a class~\code{p\_laplacian}.
501
This class interface is definied in the file \file{p\_laplacian.h}
501
This class interface is defined in the file \file{p\_laplacian.h}
502
502
and its implementation in~\file{p\_laplacian.icc}
503
503
\cindex{form!weighted!tensorial weight}%
504
504
The residual term $F(u_h)$ is computed by the member function \code{residual} while the resolution
547
547
\cindex{convergence!super-linear}
548
548
This is the major advantage of the Newton method.
549
549
Also the algorithm converge when $p\geq 3$, until $p\approx 4$.
550
\cindex{algorithm!fixed-point}%
550
\cindex{method!fixed-point}%
551
551
It was not the case with the fixed point algorithm that diverges in that case.
552
552
Finally, the Newton algorithm diverges for small values
553
553
of $p$, e.g. $p<1.5$ and the plot is not showed here.
560
560
\section{The damped Newton algorithm}
561
561
% ----------------------------------
562
562
\subsection{Principe of the algorithm}
563
\cindex{algorithm!Newton!damped}%
563
\cindex{method!Newton!damped}%
564
564
The Newton algorithm diverges when
565
565
the initial $u^{(0)}$ is too far from a solution.
566
566
Our aim is to modify the Newton algorithm and to obtain
567
a {\em globaly convergent algorith},
567
a {\em globally convergent algorithm},
568
568
i.e to converge to a solution for any initial $u^{(0)}$.
569
569
The basic idea is to decrease the step length while
570
570
maintaining the direction
629
629
&=& \langle C^{-1}F(u+\lambda \delta u) ,\
630
630
F'(u+\lambda \delta u)C^{-1}\delta u \rangle_{V,V'}
632
\cindex{operator!ajoint}%
632
\cindex{operator!adjoint}%
633
633
where the superscript $^*$ denotes the adjoint operator,
634
634
i.e. the transpose matrix the in finite dimensional case.
635
635
The practical algorithm for obtaining $\lambda$
636
636
was introduced first in~\cite{DenSch-1983} and
637
637
is also presented in~\cite[p.~385]{PreTeuVetFla-1997-C-2ed}.
638
The step length $\lambda$ that satify $\eqref{eq-backtracking}$
638
The step length $\lambda$ that satisfy $\eqref{eq-backtracking}$
639
639
is computed by using a finite sequence
640
640
$\lambda_k$, $k=0,1\ldots$ with a second order recurrence:
642
\item $k=0$~: initialisation $\lambda_0=1$.
643
If \eqref{eq-backtracking} is satified
644
whith $u + \lambda_0 \, d$ then
642
\item $k=0$~: initialization $\lambda_0=1$.
643
If \eqref{eq-backtracking} is satisfied
644
with $u + \lambda_0 \, d$ then
645
645
let $\lambda:=\lambda_0$ and the sequence stop here.
647
647
\item $k=1$~: first order recursion.
648
648
The quantities $g(0)=f(u)$ et $g'(0)=\langle f'(u),\, d\rangle$
649
are already computed at initialisation.
649
are already computed at initialization.
650
650
Also, we already have computed $g(1)=f(u+d)$ when
651
verifiying whether~\eqref{eq-backtracking} was satified.
651
verifying whether~\eqref{eq-backtracking} was satisfied.
652
652
Thus, we consider the following approximation
653
of $g(\lambda)$ by a second order polynom:
653
of $g(\lambda)$ by a second order polynomial:
655
655
\tilde{g}_1(\lambda)
661
661
After a short computation,
662
we find that the minimum of this polynom is:
662
we find that the minimum of this polynomial is:
664
664
\tilde{\lambda}_1
666
666
\Frac{-g'(0)}{2 \{ g(1) - g(0) - g'(0)\} }
668
Since the initialisation at $k=0$ does not satisfy
668
Since the initialization at $k=0$ does not satisfy
669
669
\eqref{eq-backtracking}, it is possible to show that,
670
when $\alpha$ is small enougth, we have $\tilde{\lambda}_1 \leq 1/2$
670
when $\alpha$ is small enough, we have $\tilde{\lambda}_1 \leq 1/2$
671
671
and $\tilde{\lambda}_1 \approx 1/2$.
672
672
Let $ \lambda_1 := \max(\lambda_{\rm min},\tilde{\lambda}_1)$.
673
673
If \eqref{eq-backtracking} is satisfied
677
677
\item $k\geq 2$~: second order recurrence.
678
678
The quantities $g(0)=f(u)$ et $g'(0)=\rangle f'(u),\, d\langle$
679
are available, ytogether with
679
are available, together with
680
680
$\lambda_{k-1}$, $g(\lambda_{k-1})$,
681
681
$\lambda_{k-2}$ and $g(\lambda_{k-2})$.
682
Then, $g(\lambda)$ is approximed by the following third order
682
Then, $g(\lambda)$ is approximated by the following third order
685
685
\tilde{g}_{k}(\lambda)
727
727
The sequence $(\lambda_k)_{k\geq 0}$ is strictly
728
decreasing: when the stopping criteria is not satified
728
decreasing: when the stopping criteria is not satisfied
729
729
until $\lambda_k$ reaches the machine precision $\varepsilon_{\rm mach}$
730
730
then the algorithm stops with an error.
772
772
\includegraphics{p-laplacian-damped-newton-square-r1.pdf}
775
\caption{The damped Newton algorithm on the $p$-laplacian for $d=2$:
775
\caption{The damped Newton algorithm on the $p$-Laplacian for $d=2$:
776
776
(a) convergence when $p<2$;
778
778
\label{fig-p-laplacian-damped-newton-rate}
833
833
%% \section{The affine-invariant damped Newton algorithm}
834
834
%% % ----------------------------------------------------
835
835
%% \subsection{Principe of the algorithm}
836
%% \cindex{algorithm!affine-invariant damped Newton}%
836
%% \cindex{method!affine-invariant damped Newton}%
838
838
%% The so-called {\em natural monotonicity criterion}~\cite{Deu-2004}
839
839
%% corresponds to the choice of the nonlinear preconditioner
added
removed
doc/pusrman/p-laplacian.tex
