# Solving lagrange multiplier problem via iterative solver **URL:** \ **Category:** General\ **Created:** [July 7, 2025, 5:54pm UTC](https://fenicsproject.discourse.group/t/solving-lagrange-multiplier-problem-via-iterative-solver/17914 "2025-07-07T17:54:32Z")\ **Posts on this page:** 2\ **Page:** 1
**Author:** ![kalosu](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/kalosu/32/7000_2.png) [@kalosu](https://fenicsproject.discourse.group/u/kalosu)\ **Post date:** [July 7, 2025, 5:54pm UTC](https://fenicsproject.discourse.group/t/solving-lagrange-multiplier-problem-via-iterative-solver/17914/1 "2025-07-07T17:54:32Z")
Hello there friends from the fenicsx community, I am following the example given in: [scifem/examples/real\_function\_space.py at main · scientificcomputing/scifem · GitHub](https://github.com/scientificcomputing/scifem/blob/main/examples/real_function_space.py) to solve a Neumann-only boundary value problem by means of a Lagrangian multiplier. In general, if I follow the example provided in the link, the results that I obtain work nicely. However, once I try to scale my simulation to a denser grid, the steps related to solving this linear system take a significant amount of my computation time. I was wondering if it was possible to use an iterative solver in combination with this realfunction space approach to speed up the evaluation process? As mentioned above, I follow exactly the provided example, even with the settings for the direct solver. Any comment will be highly appreciated. ---
**Author:** ![dparsons](https://avatars.discourse-cdn.com/v4/letter/d/5e9695/32.png) [@dparsons](https://fenicsproject.discourse.group/u/dparsons)\ **Post date:** [July 9, 2025, 9:23am UTC](https://fenicsproject.discourse.group/t/solving-lagrange-multiplier-problem-via-iterative-solver/17914/2 "2025-07-09T09:23:13Z")
A request for this sort of capability was recently requested on the scifem repo at [Real elements with non-linear problems · Issue #114 · scientificcomputing/scifem · GitHub](https://github.com/scientificcomputing/scifem/issues/114) but the requestor found that a possible solution was already available at > > > Suggestion by @CastillonMiguel (and worked out with him) > > \`\`\`python > \# # Real fun…ction spaces > \# > \# Author: Jørgen S. Dokken > \# > \# License: MIT > \# > \# In this example we will show how to use the "real" function space to solve > \# a singular Poisson problem. > \# The problem at hand is: > \# Find $u \\in H^1(\\Omega)$ such that > \# \\begin{align} > \# -\\Delta u &= f \\quad \\text{in } \\Omega, \\\\ > \# \\frac{\\partial u}{\\partial n} &= g \\quad \\text{on } \\partial \\Omega, \\\\ > \# \\int\_\\Omega u &= h. > \# \\end{align} > \# > \# ## Lagrange multiplier > \# We start by considering the equivalent optimization problem: > \# Find $u \\in H^1(\\Omega)$ such that > \# \\begin{align} > \# \\min\_{u \\in H^1(\\Omega)} J(u) = \\min\_{u \\in H^1(\\Omega)} \\frac{1}{2}\\int\_\\Omega \\vert \\nabla u \\cdot \\nabla u \\vert \\mathrm{d}x - \\int\_\\Omega f u \\mathrm{d}x - \\int\_{\\partial \\Omega} g u \\mathrm{d}s, > \# \\end{align} > \# such that > \# \\begin{align} > \# \\int\_\\Omega u = h. > \# \\end{align} > \# We introduce a Lagrange multiplier $\\lambda$ to enforce the constraint: > \# \\begin{align} > \# \\min\_{u \\in H^1(\\Omega), \\lambda\\in \\mathbb{R}} \\mathcal{L}(u, \\lambda) = \\min\_{u \\in H^1(\\Omega), \\lambda\\in \\mathbb{R}} J(u) + \\lambda (\\int\_\\Omega u \\mathrm{d}x-h). > \# \\end{align} > \# We then compute the optimality conditions for the problem above > \# \\begin{align} > \# \\frac{\\partial \\mathcal{L}}{\\partial u}\[\\delta u\] &= \\int\_\\Omega \\nabla u \\cdot \\nabla \\delta u \\mathrm{d}x + \\lambda\\int \\delta u \\mathrm{d}x - \\int\_\\Omega f \\delta u ~\\mathrm{d}x - \\int\_{\\partial \\Omega} g \\delta u~\\mathrm{d}s = 0, \\\\ > \# \\frac{\\partial \\mathcal{L}}{\\partial \\lambda}\[\\delta \\lambda\] &=\\delta \\lambda (\\int\_\\Omega u \\mathrm{d}x -h)= 0. > \# \\end{align} > \# We write the weak formulation: > \# > \# $$ > \# \\begin{align} > \# \\int\_\\Omega \\nabla u \\cdot \\nabla \\delta u~\\mathrm{d}x + \\int\_\\Omega \\lambda \\delta u~\\mathrm{d}x = \\int\_\\Omega f \\delta u~\\mathrm{d}x + \\int\_{\\partial \\Omega} g v \\mathrm{d}s\\\\ > \# \\int\_\\Omega u \\delta \\lambda \\mathrm{d}x = h \\delta \\lambda . > \# \\end{align} > \# $$ > \# > \# where we have moved $\\delta\\lambda$ into the integral as it is a spatial constant. > > \# ## Implementation > \# We start by creating the domain and derive the source terms $f$, $g$ and $h$ from our manufactured solution > \# For this example we will use the following exact solution > \# \\begin{align} > \# u\_{exact}(x, y) = 0.3y^2 + \\sin(2\\pi x). > \# \\end{align} > > from mpi4py import MPI > from petsc4py import PETSc > import dolfinx.fem.petsc > > import numpy as np > from scifem import create\_real\_functionspace, assemble\_scalar > import ufl > > M = 20 > mesh = dolfinx.mesh.create\_unit\_square( > MPI.COMM\_WORLD, M, M, dolfinx.mesh.CellType.triangle, dtype=np.float64 > ) > V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1)) > > > def u\_exact(x): > return 0.3 \* x\[1\] \*\* 2 + ufl.sin(2 \* ufl.pi \* x\[0\]) > > > x = ufl.SpatialCoordinate(mesh) > n = ufl.FacetNormal(mesh) > g = ufl.dot(ufl.grad(u\_exact(x)), n) > f = -ufl.div(ufl.grad(u\_exact(x))) > h = assemble\_scalar(u\_exact(x) \* ufl.dx) > > \# We then create the Lagrange multiplier space > > R = create\_real\_functionspace(mesh) > > \# Next, we can create a mixed-function space for our problem > > > W = ufl.MixedFunctionSpace(V, R) > u = dolfinx.fem.Function(V) > lmbda = dolfinx.fem.Function(R) > w = \[u, lmbda\] > du, dl = ufl.TestFunctions(W) > dw = ufl.TrialFunctions(W) > \# We can now define the variational problem > > zero = dolfinx.fem.Constant(mesh, dolfinx.default\_scalar\_type(0.0)) > > F = ufl.inner(ufl.grad(u), ufl.grad(du)) \* ufl.dx > F += ufl.inner(lmbda, du) \* ufl.dx > F += ufl.inner(u, dl) \* ufl.dx > F -= ufl.inner(f, du) \* ufl.dx + ufl.inner(g, du) \* ufl.ds > F += ufl.inner(zero, dl) \* ufl.dx > F\_blocked = ufl.extract\_blocks(F) > > > J = sum(ufl.derivative(F, w\[i\], dw\[i\]) for i in range(len(w))) > J\_blocked = ufl.extract\_blocks(J) > > import scifem > class RealSpaceNewtonSolver(scifem.BlockedNewtonSolver): > > def \_assemble\_residual(self, x: PETSc.Vec, b: PETSc.Vec) -\> None: > """Assemble the residual F into the vector b. > Args: > x: The vector containing the latest solution > b: Vector to assemble the residual into > """ > # Assemble F(u\_{i-1}) - J(u\_D - u\_{i-1}) and set du|\_bc= u\_D - u\_{i-1} > with b.localForm() as b\_local: > b\_local.set(0.0) > dolfinx.fem.petsc.assemble\_vector\_block( > b, > self.\_F, > self.\_a, > bcs=self.\_bcs, > x0=x, > alpha=-1.0, > ) > start\_pos = 0 > for s in self.\_u: > num\_sub\_dofs = ( > s.function\_space.dofmap.index\_map.size\_local > \* s.function\_space.dofmap.index\_map\_bs > ) > # FIXME: Add documentation here! > if s.function\_space.dofmap.index\_map.size\_global == 1: > assert s.function\_space.dofmap.index\_map\_bs == 1 > if s.function\_space.dofmap.index\_map.size\_local \> 0: > b.array\_w\[start\_pos:start\_pos+num\_sub\_dofs\] -= h > start\_pos += num\_sub\_dofs > b.ghostUpdate(PETSc.InsertMode.INSERT\_VALUES, PETSc.ScatterMode.FORWARD) > > \# Note that we have defined the variational form in a block form, and > \# that we have not included $h$ in the variational form. We will enforce this > \# once we have assembled the right hand side vector. > > solver = RealSpaceNewtonSolver(F\_blocked, w,J=J\_blocked, bcs=\[\], > petsc\_options={"ksp\_type": "preonly", "pc\_type": "lu", "pc\_factor\_mat\_solver\_type": "mumps"}) > solver.solve() > > print(f"{dolfinx.fem.assemble\_scalar(dolfinx.fem.form(w\[0\]\*ufl.dx))}") > diff = w\[0\] - u\_exact(x) > error = dolfinx.fem.form(ufl.inner(diff, diff) \* ufl.dx) > > print(f"L2 error: {np.sqrt(assemble\_scalar(error)):.2e}") > > with dolfinx.io.XDMFFile(MPI.COMM\_WORLD, "solution.xdmf", "w") as xdmf: > xdmf.write\_mesh(mesh) > xdmf.write\_function(w\[0\]) > > \`\`\`