p-Laplace equation

Potential for Laplace equation

Task 1

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

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!

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[0])*sin(2.*pi*x[1])", 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[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.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[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.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()