Tutorial 25: Lagrange multipliers

In this tutorial, we will learn

How to enforce constraints using Lagrange multipliers
How to work with ConstantFESpace

Problem statement

In this tutorial, we solve the Poisson equation with pure Neumann boundary conditions. This problem is well-known to be singular since the solution is defined up to a constant, which we have to fix to obtain a unique solution. Here, we will use a Lagrange multiplier to enforce that the mean value of the solution equals a given constant.

The problem reads: find $u$ and $λ$ such that

\[\left\lbrace \begin{aligned} -\Delta u = f \ &\text{in} \ \Omega,\\ \nabla u\cdot n = g \ &\text{on}\ \Gamma,\\ \int_{\Omega} u \ {\rm d}\Omega = \bar{u},\\ \end{aligned} \right.\]

where $\Omega$ is our domain, $\Gamma$ is its boundary, $n$ is the outward unit normal vector, and $\bar{u}$ is a given constant that fixes the mean value of the solution.

Numerical scheme

The weak form of this problem using Lagrange multipliers reads: find $(u,λ) \in V \times \Lambda$ such that

\[\begin{aligned} \int_{\Omega} \nabla u \cdot \nabla v \ {\rm d}\Omega + \int_{\Omega} λv \ {\rm d}\Omega + \int_{\Omega} uμ \ {\rm d}\Omega = \int_{\Omega} fv \ {\rm d}\Omega + \int_{\Gamma} v(g\cdot n) \ {\rm d}\Gamma + \int_{\Omega} μ\bar{u} \ {\rm d}\Omega \end{aligned}\]

for all $(v,μ) \in V \times \Lambda$, where $V = H^1(\Omega)$ and $\Lambda = \mathbb{R}$.

Implementation

First, we load the Gridap package and define the exact solution that we will use to manufacture the source term and boundary condition:

using Gridap

u_exact(x) = sin(x[1]) * cos(x[2])

Now we can create a simple Cartesian mesh of the unit square:

model = CartesianDiscreteModel((0,1,0,1),(8,8))

We will use first order Lagrangian finite elements for the primal variable u.

order = 1
reffe = ReferenceFE(lagrangian, Float64, order)
V = FESpace(model, reffe)

For the Lagrange multiplier λ, we need a space of constant functions, since λ ∈ ℝ. In Gridap, we can create such a space using ConstantFESpace:

Λ = ConstantFESpace(model)

Conceptually, a ConstantFESpace is a space defined on the whole domain with a single degree of freedom, which is what we need for the Lagrange multiplier λ. We finally bundle both spaces into a multi-field space:

X = MultiFieldFESpace([V, Λ])

Integration

We need to create the triangulation and measures for both domain and boundary integration:

Ω = Triangulation(model)
Γ = BoundaryTriangulation(model)
dΩ = Measure(Ω, 2*order)
dΓ = Measure(Γ, 2*order)

Next, we manufacture the source term f and Neumann boundary condition g from the exact solution. We also compute the mean value ū that we want to enforce:

f(x) = -Δ(u_exact)(x)
g(x) = ∇(u_exact)(x)
ū = sum(∫(u_exact)dΩ)
nΓ = get_normal_vector(Γ)

Weak Form

We can now define the bilinear and linear forms of our problem. Note how the forms take tuples as arguments, representing the multi-field nature of our solution:

a((u,λ),(v,μ)) = ∫(∇(u)⋅∇(v) + λ*v + u*μ)dΩ
l((v,μ)) = ∫(f*v + μ*ū)dΩ + ∫(v*(g⋅nΓ))*dΓ

Solution

We can now create the FE operator and solve the system:

op = AffineFEOperator(a, l, X, X)
uh, λh = solve(op)

Note how we get two values from solve: the primal solution uh and the Lagrange multiplier λh. Finally, we compute the L2 error and verify that the mean value constraint is satisfied:

eh = uh - u_exact
l2_error = sqrt(sum(∫(eh⋅eh)*dΩ))
ūh = sum(∫(uh)*dΩ)

The L2 error should be small (of order h²) and ūh should be very close to ū, showing that both the equation and the constraint are well satisfied.

Visualization

We can visualize the solution and error by writing them to a VTK file:

writevtk(Ω, "results", cellfields=["uh"=>uh, "error"=>eh])

This page was generated using Literate.jl.