Variational formulation using Lagrange multiplier

variational formulation
1

I tried to understand the implementation of lagrange multipler through https://fenicsproject.org/olddocs/dolfin/latest/python/demos/neumann-poisson/demo_neumann-poisson.py.html

I don’t understand how the variational form takes c and d terms in bilinear form. Can you give some explanation of how you formed this equation. This seems to be different from general lagrange multiplier applications.

My major purpose is to solve the periodic homogenization problem using Lagrange multiplier.
https://comet-fenics.readthedocs.io/en/latest/demo/periodic_homog_elas/periodic_homog_elas.html

2

If you’re familiar with Lagrange multipliers in the context of constrained minimization, the easiest way to derive the variational form from the linked example is by starting with the following energy functional for the Poisson equation with Neumann boundary conditions:

\mathcal{J}(u) = \frac{1}{2}\int_\Omega\vert\nabla u\vert^2\,d\Omega - \int_\Omega fu\,d\Omega - \int_{\Gamma_N}gu\,d\partial\Omega\text{ .}

If the Neumann boundary part \Gamma_N were only part of \partial\Omega, and there were another part where u was constrained to be zero, then the Poisson problem could be solved uniquely by minimizing this energy functional. However, when \Gamma_N = \partial\Omega, the solution is only defined up to a constant; under the compatibility condition

\int_\Omega f\,d\Omega = -\int_{\partial\Omega}g\,d\partial\Omega\text{ ,}

we can add a constant to u without changing \mathcal{J}(u). Thus, we want to add the constraint

\int_\Omega u\,d\Omega = 0\text{ ,}

to pin down this free constant. To do this using the method of Lagrange multipliers, you would seek a stationary point of the Lagrangian

\mathcal{L}(u,\lambda) = \mathcal{J}(u) + \lambda\int_\Omega u\,d\Omega\text{ .}

Taking a variation (i.e., the Gateaux derivative in the direction of an arbitrary function/multiplier pair (\delta u,\delta\lambda)) and setting it to zero, we get

\delta\mathcal{L}(u,\lambda) = \int_\Omega\nabla u\cdot\nabla\delta u\,d\Omega - \int_\Omega f\delta u\,d\Omega - \int_{\Gamma_N}g\delta u\,d\partial\Omega + \delta\lambda\int_\Omega u\,d\Omega + \lambda\int_\Omega\delta u\,d\Omega = 0\text{ .}

If \lambda and \delta\lambda are in \mathbb{R}, we can pull them inside the integrals over \Omega. Renaming \delta u\to v, \lambda\to c, and \delta\lambda \to d, we recover the formulation from the linked tutorial.

(I’ll add a brief note that you can also use multipliers to apply constraints to weak problems that do not correspond to minimization, by just directly adding the variations of the multiplier terms to a weak form (e.g., the pressure in the weak form of the incompressible Navier–Stokes equations, which acts as a Lagrange multiplier for the incompressibility constraint). But, for whatever reason, Lagrange multipliers are most frequently taught in the context of minimization problems.)

4 Likes
3

I do not understand why you choose that particular form of constrain for u. Can you elaborate why you choose that particular constrain? Could I choose another constrain? If so, can you give me another example of another constrain I can choose?

Thanks,
Jorge

4

@jormoral, The constrain depend on your application. For above case, it is one of many options to make u particular to constrain the rigid body motion (for above, it might be translation). You can fix u at any specific node/location and that can also serve as one option for constraining.

1 Like
5

I am trying to understand how the above approach translates to a case where one wants to constrain the flux for a set of problems such as \nabla\cdot(\kappa \nabla u)=0 in \Omega_1\bigcup\Omega_2 and \nabla\cdot(D\nabla c)=0 in \Omega_2 where the constraint at the dofs at \Gamma_X = \partial\Omega_1\bigcup\partial\Omega_2 is \lambda\Bigl(\kappa\nabla u\cdot \pmb{n} - FD\nabla c\cdot \pmb{n}\Bigr)=0.

I’m unclear how the Lagrange\cdot \pmb{n} constraint gets added to the energy functional since it is not integrated like the demo.

\mathcal{J}(u, c) = \frac{1}{2}\int_{\Omega}\vert \nabla u \vert^2\mathrm{dx} - \int_{\Gamma_X}g u\mathrm{ds} + \frac{1}{2}\int_{\Omega_2}D\vert \nabla c \vert^2\mathrm{dx} - \int_{\Gamma_X}g_2c\mathrm{ds}
6

We add Lagrange multiplier term to energy as like adding zero. If you see the coefficient of \lambda (given constraint) is actually zero which doesn’t affect the energy. However, the benefit of adding is that we add extra dofs (\lambda), to make the constrained previous problem with boundary conditions to unconstrained problem. If there are 4 constraints applied, the new dof would be (ndof+4).

Unconstrained problem has characteristic that, the linear form Ax=b with A_{(ndof+4,ndof+4)} and b_{ndof+4} has constrained already involved in it, and can be directly solved to get the required x. You can find more detail about lagrange in calculus of variations references.

1 Like
7

Thank you. Because the variation form directly from the Lagrangian introduced zeros on diagonals I opted for a perturbed Lagrangian:

\delta\mathcal{L} = \int_{\Omega}\kappa\nabla u \nabla \delta u\mathrm{dx} - \int_{\partial\Omega}g\delta u\mathrm{ds} + \int_{\Omega_2}D\nabla c \nabla \delta c\mathrm{dx} - \int_{\partial\Omega_1\bigcup\partial\Omega_2}g_2\delta c\mathrm{ds} + \int_{\partial\Omega_1\bigcup\partial\Omega_2}\delta\lambda\bigl(\kappa\nabla u\cdot\pmb{n} - FD\nabla c \cdot \pmb{n}\bigr)\mathrm{ds} + \int_{\partial\Omega_1\bigcup\partial\Omega_2}\frac{1}{\alpha}\lambda\delta\lambda\mathrm{ds}

for sufficiently large \alpha\in\mathbb{R}^+.

However, I could use dolfiny restriction GitHub - michalhabera/dolfiny: Dolfin-y, high level wrappers for dolfin-x, the FEniCS library to get rid of the zero diagonals.

I would like to extend the problem to solve heat equation in \Omega_2 coupled to Poisson equation in \Omega=\Omega_1\bigcup\Omega_2 but adding \int_{\Omega_2}\frac{c}{\Delta t} \delta c \mathrm{dx} - \int_{\Omega_2}\frac{c_0}{\Delta t} \delta c \mathrm{dx} for initial condition c(t=0)=c_0 and time step \Delta t does not seem to be correct.

\delta\mathcal{L} = \int_{\Omega}\kappa\nabla u \nabla \delta u\mathrm{dx} - \int_{\partial\Omega}g\delta u\mathrm{ds} + \int_{\Omega_2}D\nabla c \nabla \delta c\mathrm{dx} - \int_{\partial\Omega_1\bigcup\partial\Omega_2}g_2\delta c\mathrm{ds} + \int_{\partial\Omega_1\bigcup\partial\Omega_2}\delta\lambda\bigl(\kappa\nabla u\cdot\pmb{n} - FD\nabla c \cdot \pmb{n}\bigr)\mathrm{ds} + \int_{\partial\Omega_1\bigcup\partial\Omega_2}\frac{1}{\alpha}\lambda\delta\lambda\mathrm{ds} + \int_{\Omega_2}\frac{c}{\Delta t} \delta c \mathrm{dx} - \int_{\Omega_2}\frac{c_0}{\Delta t} \delta c \mathrm{dx}