Choosing optimization procedure for solving non-linear PDE with finite-element method

17 Views Asked by At

I'm trying to write a code solving the Dirichlet boundary problem for $p$-Laplacian on an arbitrary planar domain $\Omega$:

$$\begin{cases} -\Delta_p u = f \text{ on } \Omega, \\ u\big|_{\partial \Omega} \equiv 0 \end{cases} $$

using finite-elements method. I rely on the classical formulation as in Barret, Liu 1993: https://www.ams.org/journals/mcom/1993-61-204/S0025-5718-1993-1192966-4/S0025-5718-1993-1192966-4.pdf: given a finite-element continuous piecewise-linear basis functions, minimize

$$ J(u) = \frac{1}{p} \int_{\Omega} |\nabla u|^p \: d\Omega - \int_{\Omega} fv \: d\Omega $$

For the optimization routine I decided to implement some form of Polak-Rebiere conjugate gradient method with global convergence (since PR is known to have been already used for this particular problem), and tried programming this one: http://dx.doi.org/10.1080/01630568908816355 It mostly works on typical test functions, although sometimes it doesn't converge to the global minimum if the initial value is too far from the minimum (in such cases it's usually possible to make the method converge by some changes on constant parameters it uses). It doesn't seem to work for my functional optimization though (and I don't know if it is the gradient method or something else, for example integrals computation, causing the problems). As I understand, the global convergence of such approaches can be guaranteed only if the initial approximation $x_0$ is chosen correctly (for example, the set $\{ x: J(x) \leqslant x_0 \}$ must be at least bounded), so I wonder if there's any way to assure that the selected optimization method will converge from any initial value $x_0$, or to choose an initial approximation to guarantee convergence? What is a standart way of doing that? Would be greatful for any advice

P.S. Here's a scheme for the algorithm I've tried to implement. Maybe there're more reliable forms of conjugate gradient method to use? Do I get it correctly, that $\alpha_k$ could be computed by finding the least power $j\in \{0, 1, \dots\}$ of $t$, for which $\alpha_k = \rho_k t^j$ satisfies (22-23)?