# p-Laplace equation¶

## Potential for Laplace equation¶

Formulate Laplace equation

\begin{align}\begin{aligned}-\Delta u &= f &&\text{ in } \Omega,\\ u &= 0 &&\text{ on } \partial\Omega,\end{aligned}\end{align}

as a 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¶

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$$.

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$$.

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!

## Reference solution¶

Reference solution consists of three files – one module file:

Show/Hide File p_laplace.py

Download Code

from dolfin import *
import matplotlib.pyplot as plt

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

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)

plt.gcf().show()
plt.clf()
plot(u, mode="warp", title='p-Laplace, p=%g, eps=%g'%(float(p), float(eps)))
plt.gcf().canvas.draw()

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

return u, energy_regularized, energy


and two executable scripts:

Show/Hide File p_small.py

Download Code

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

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
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.figure()
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()


Show/Hide File p_large.py

Download Code

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

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
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.figure()
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()