p-Laplace equation

Potential for Laplace equation

Task 1. Formulate Laplace equation

\[\begin{split}-\Delta u &= f \quad \text{ in } \Omega, \\ u &= 0 \quad \text{ on } \partial\Omega,\end{split}\]

as variational problem (i.e. find potential for the equation) and solve it using FEniCS with data

• \(\Omega = [0, 1] \times [0, 1]\),
• \(f = 1 + \cos(2 \pi x) \sin(2 \pi y)\).

Use UFL function derivative to automatically obtain Galerkin formulation from the potential. Don’t assume linearity of the PDE - solve it as nonliner (Newton will converge in 1 step).

Potential for p-Laplace equation

Task 2. Replace every occurrence of number \(2\) in potential for Laplace equation by \(p\). This is called \(p\)-Laplacian for \(1 < p < +\infty\).

Convince yourself that resulting PDE is non-linear whenever \(p \neq 2\).

Task 3. Run the algorithm from Task 1 with \(p = 1.1\) and \(p = 11\).

Hint. Use DOLFIN class Constant to avoid form recompilation by FFC for distinct values of \(p\).

Do you know where is the problem? If not, compute second Gateaux derivative of the potential (which serves as Jacobian for the Newton-Rhapson algorithm) and look at its value for \(u = 0\).

Task 4. Add some regularization to the potential to make it non-singular/non-degenerate. (Prefer regularization without square roots.) Find a solution of the original \(p\)-Laplace problem with \(p = 1.1\) and \(p = 11\) using careful approximation by regularized problem. Report \(max_\Omega u_h\) for \(u_h\) being the approximate solution.

Hint. For u_h Lagrange degree 1 Function the maximum matches maximal nodal value so it is u_h.vector().max() because nodal basis is used. Warning. This does not generally hold for higher polynomial degrees!

Note

Lecturer note. Totally wrong solution may be easily obtained by overregularization. Student must find a (heuristic) way to show that convergence is established.

The trick in regularization-convergence for \(p = 11\) is that Newton solver needs to be initialized by previous result computed using higher regularization - see for-loop in p_large.py.

Convergence in discretization parameter is not so important here compared to regularization and good solution is obtained on mild meshes. Also exact, appropriately high quadrature for large p is not needed here - parameters['form_compiler']['quadrature_degree'] = 11 yields the same results.

Solution for \(p = 1.1,\,11\) has maximum around 1.3E-7 and 0.41 respectively.

Reference solution

File p_laplace.py

from dolfin import *

mesh = UnitSquareMesh(40, 40)
V = FunctionSpace(mesh, 'Lagrange', 1)
f = Expression("1.+cos(2.*pi*x[0])*sin(2.*pi*x[1])")

def p_laplace(p, eps, u0=None):
    """Solves regularized p-Laplacian with mesh, space and right-hand side
    defined above. Returns solution, regularized energy and energy."""
    p = Constant(p)
    eps = Constant(eps)

    # Initial approximation for Newton
    u = u0.copy(deepcopy=True) if u0 else Function(V)

    # Energy functional
    E = ( 1./p*(eps + dot(grad(u), grad(u)))**(0.5*p) - f*u ) * dx

    # First Gateaux derivative
    F = derivative(E, u)

    # Solve nonlinear problem; Newton is used by default
    bc = DirichletBC(V, 0.0, lambda x,onb: onb)
    solver_parameters = {'newton_solver': {'maximum_iterations': 1000}}
    solve(F == 0, u, bc, solver_parameters=solver_parameters)

    plot(u, title='p-Laplace, p=%g, eps=%g'%(float(p), float(eps)))

    # Compute energies
    energy_regularized = assemble(E)
    eps.assign(0.0)
    energy = assemble(E)

    return u, energy_regularized, energy

File p_small.py

from p_laplace import p_laplace
import numpy as np
import matplotlib.pyplot as plt
from dolfin import interactive

epsilons = [10.0**i for i in range(-10, -22, -2)]
energies = []
maximums = []

# Converge with regularization
for eps in epsilons:
    result = p_laplace(1.1, eps)
    u = result[0]
    energies.append(result[1:])
    # Maxmimal nodal value (correct maximum only for P1 function)
    maximums.append(u.vector().max())
energies = np.array(energies)

# Report
print ' epsilon  |  energy_reg  |  energy | u max \n', \
    np.concatenate( ( np.array([epsilons,]).T,
                      energies,
                      np.array([maximums,]).T,
                    ), axis=1)

# Plot energies
plt.plot(epsilons, energies[:,0], 'o-', label='energy regularized')
plt.plot(epsilons, energies[:,1], 'o-', label='energy')
plt.xscale('log')
plt.legend(loc='upper left')
plt.show(block=False)

interactive()

File p_large.py

from p_laplace import p_laplace
import numpy as np
import matplotlib.pyplot as plt
from dolfin import interactive

epsilons = [10.0**i for i in np.arange(1.0, -6.0, -0.5)]
energies = []
maximums = []

# Converge with regularization,
# use previous solution to start next Newton!
u = None
for eps in epsilons:
    result = p_laplace(11.0, eps, u)
    u = result[0]
    energies.append(result[1:])
    # Maxmimal nodal value (correct maximum only for P1 function)
    maximums.append(u.vector().max())
energies = np.array(energies)

# Report
print ' epsilon  |  energy_reg  |  energy | u max \n', \
    np.concatenate( ( np.array([epsilons,]).T,
                      energies,
                      np.array([maximums,]).T,
                    ), axis=1)

# Plot energies
plt.plot(epsilons, energies[:,0], 'o-', label='energy regularized')
plt.plot(epsilons, energies[:,1], 'o-', label='energy')
plt.xscale('log')
plt.legend(loc='upper left')
plt.show(block=False)

interactive()