# Poisson equation: Impose full gradient as boundary condition via Lagrange multipliers

I have a physical problem governed by the Poisson equation in two dimensions $$-\nabla^2 u = f(x,y), \; in \; \Omega$$ I have measurements of the two gradient components $\partial{u}/\partial{x}$ and $\partial{u}/\partial{y}$ along some part of the boundary, $\Gamma_m$, so would like to impose $$\frac{\partial{u}}{\partial{x_i}}_0 = g_m, \; on \; \Gamma_m$$ and propagate into the far field.

The tangential gradient component, $\frac{\partial{u}}{\partial{x}}_{(t,0)}$, I can just integrate and then enforce through a Dirichlet condition, such that $$\int_{\Gamma_m}\frac{\partial{u}}{\partial{x}}_{(t,0)} \, ds = u_0$$ In order to simultaneously impose the normal component, $\frac{\partial{u}}{\partial{x}}_{(n,0)}$, I gathered I would have to go via Lagrange multipliers.

So I think the variational form is then $$\int_\Omega \nabla{u} \cdot\nabla{v}\, dx - \lambda \int_{\Gamma_m} ( \frac{\partial{u}}{\partial{x}}_{(n,0)}-g_m ) v \,ds = \int_\Omega f\, v\, dx$$ I spent a long time trying to piece it together from the information on related problems such as https://answers.launchpad.net/fenics/+question/212434 https://answers.launchpad.net/fenics/+question/216323

but still cannot see where I am going wrong. My solution attempt so far is:

from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(64, 64)
V = FunctionSpace(mesh, "Lagrange", 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R

# Create mesh function over cell facets
boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)

# Mark left boundary facets as subdomain 0
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] < DOLFIN_EPS

Gamma_Left = LeftBoundary()
Gamma_Left.mark(boundary_parts, 0)

class FarField(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and ( (x[0] > 1.0-DOLFIN_EPS) \
or (x[1]<DOLFIN_EPS) or (x[1]> 1.0-DOLFIN_EPS) )

Gamma_FF = FarField()
Gamma_FF.mark(boundary_parts, 1)

# Define boundary condition
u0 = Expression("sin(x[1]*pi)")
bcs = [DirichletBC(V, u0, Gamma_Left)]

# Define variational problem
(u, lmbd) = TrialFunctions(W)
(v, d) = TestFunctions(W)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Constant(0.0)
h = Constant(-4.0)
n = FacetNormal(mesh)

(f*v*dx + g*v*ds(1) + h*d*ds(0) + lmbd*h*ds(0))

a = lhs(F)
L = rhs(F)

# Compute solution
A = assemble(a, exterior_facet_domains=boundary_parts)
b = assemble(L, exterior_facet_domains=boundary_parts)
for bc in bcs: bc.apply(A, b)

w = Function(W)
solve(A, w.vector(), b, 'lu')
(u,lmbd) = w.split()

# Plot solution
plot(u, interactive=True)


which runs but gives a noisy result not at all resembling a solution to a Poisson equation. It seems to have something to do with the combined function spaces, but I cannot find the mistake.
I would appreciate any help or pointers in the right direction - many thanks already!
Cheers
Markus

• Let me get this right: You have both Dirichlet and Neumann data, but only on part of the boundary? – Christian Clason May 15 '13 at 14:46
• As I have understood the OP, it is the gradient that is given at the boundary. The Dirichlet data is used to impose the tangential derivative. I thought it's odd to impose both Dirichlet and Neumann at one part of the boundary but maybe in this particular situation it is consistent. So, the problem is rather how to apply gradient data at the boundary (via multipliers). – Jan May 15 '13 at 21:25
• True, that would give consistent data, but you still have the problem of lack of stability, and the fact that you have boundary conditions on only part of the boundary. – Christian Clason May 16 '13 at 6:19
• Ok, let me give some more info on the specific physical problem I am trying to solve. I have a static magnetic field that I can reasonably assume to be rotationally symmetric, thus 2D. I measure radial and axial components of magnetic field density vector along a line, fairly close to the axis of rotation and would like to see this magnetic field at substantial distance away from this axis of rotation. The combination of Dirichlet and Neumann BC was just my idea of approaching the problem as Jan eloquently described - imposed gradient data at the boundary. – Markus May 16 '13 at 7:15
• OK, that changes things significantly. So you have an unbounded domain, and derivative information on the whole "finite" part of the boundary? – Christian Clason May 16 '13 at 7:35

First, a general point: you cannot prescribe arbitrary boundary conditions for a partial differential operator and expect that the partial differential equation (which always includes both operator and boundary conditions) is well-posed, i.e., admits a unique solution that depends continuously on the data -- all of which is a necessary condition for actually trying to compute something.

Depending on the operator, there are often quite a number of valid conditions you can impose (to get a taste, see the three-volume monograph by Lions and Magenes). However, what you are trying to do (specify the full gradient, which is equivalent to both Dirichlet and Neumann conditions on the same (part of the) boundary for a second-order elliptic PDE) is not among them - this is known as a lateral Cauchy problem, and is ill-posed: there is no guarantee that a given pair of boundary data admits a solution, and even if one exists, there is no stability with respect to small perturbations in the data. (In fact, this is the original ill-posed problem in the sense of Hadamard, and the classical example why you cannot ignore boundary conditions when discussing well-posedness. You can find an explicit example in his Lectures on Cauchy's problem in linear partial differential equations from the 1920s.)

Now to your specific problem (which might be a textbook example of an XY problem). If I understand correctly, you wish to solve an exterior problem for the Poisson equation, and you take a two-dimensional computional domain $(r,R)\times(a,b)$. You have full data (Dirichlet - using the tangential derivative trick - and Neumann) on (say) $x=r$. If $R$ is sufficiently large, you can justify a radiation condition (which prescribes the asymptotic behavior of your solution a t $x\to \infty$); see, e.g., the book Numerical methods for exterior problems by Long'an Ying. The question is what you know about the remaining boundary parts, $y=a$ and $y=b$.

1. If you can impose boundary conditions (Neumann, Robin, Dirichlet - which you would need to fix the constant in the integration of the tangential derivative, by the way), then it is enough to use either the normal components of your gradient as Neumann condition (if you can fix the constant mode) or integrate the tangential components as a Dirichlet condition. Since both conditions presumably correspond to the same function, you get the same solution either way.

2. If you don't know the behavior at $y=a$ and $y=b$, you really have a lateral Cauchy problem, and you cannot compute a solution using standard finite element methods. The standard way of dealing with this is the method of quasi-reversibility (introduced by Lattés and Lions in the 1960s; see their book), which consists in approximating the second-order problem by a fourth-order problem (where you can - and need to - prescribe two boundary conditions). In your case, this would amount to replacing $-\Delta u=f$ by $-\Delta u + \varepsilon \Delta^2 u=f$ for some small $\varepsilon>0$. (This can also be interpreted as minimizing the residual in a suitable norm and adding an $H^2$ regularization term; related to your comment on treating the problem as an inverse problem.) You can show that for compatible data (i.e., a pair of boundary conditions that actually correspond to a solution $u$ of the Poisson equation), the quasi-reversibility solutions $u_\varepsilon$ converge to $u$ as $\varepsilon\to0$.

Since this is now a fourth-order problem with solutions in $H^2$, the best way of solving it numerically is using a mixed finite element formulation such as the one described in the paper by Dardé, Hannukainen and Hyvönen. (There are also some slides online.) It should not be too difficult to implement this approach using FEniCS.

• For implementation by mixed elements in FEniCS see biharmonic demo. This is probably without Laplace term but I guess it can be easily added. – Jan Blechta May 15 '13 at 18:58
• Hi Christian, thanks for your suggestion! I was under the impression that the Poisson equation was benign as far as numerical stability is concerned - thanks for pointing that out. I'll read up on it as you suggested. Not sure if this changes things substantially, but as mentioned in the comment further up, the Dirichlet-Neumann combination is perhaps misleading. 'All' I am looking for, is a way of imposing gradient data at the boundary. – Markus May 16 '13 at 7:20
• The Poisson equation is benign, but that's not the equation you're trying to solve :) (Boundary conditions are an integral part of the equation; the operator alone is insufficient to decide stability.) – Christian Clason May 16 '13 at 7:30
• Alright, that gives me something to chew on. Thanks everyone for your time, advice and patience - and my apologies for falling into the XY trap... – Markus May 20 '13 at 17:35

You can't expect that solution to your altered problem would be a solution to Poisson problem because you need to change the problem somehow to make it well-posed.

One could guess that possible problem formulation is to minimize $$F(u, \lambda) = \int_\Omega \frac{1}{2}|\nabla u|^2 \;\mathrm{d}x - \int_\Omega f u \;\mathrm{d}x - \int_{\Gamma_\mathrm{N}} g u \;\mathrm{d}S + \int_{\Gamma_\mathrm{N}} \lambda (u-u_\mathrm{D})\;\mathrm{d}S$$ over $(u,\lambda) \in V\times L^2(\Gamma_\mathrm{N})$ where $V = \{v\in H^1;v|_{\Gamma_\mathrm{D}}=0\}$ is space with zero trace on Dirichlet boundary $\Gamma_\mathrm{D}$ and $u_\mathrm{D}$ is Dirichlet value you would like to enforce additionaly on Neumann boundary $\Gamma_\mathrm{N}$. Without the last term minimization of $F(u)$ yields $$0 = \mathrm{D}F(u)[v] = \int_\Omega \nabla u\cdot\nabla v\;\mathrm{d}x - \int_\Omega fv\;\mathrm{d}x - \int_{\Gamma_\mathrm{N}} gv\;\mathrm{d}S \quad \forall v\in V,$$ which is exactly Poisson problem with Neumann condition on ${\Gamma_\mathrm{N}}$ and zero (it is simply adjusted to non-zero) Dirichlet condition on ${\Gamma_\mathrm{D}}$. With additional Lagrange multiplier term necessary condition of minimality reads $$0 = \mathrm{D}F(u,\lambda)[v,\mu] = \int_\Omega \nabla u\cdot\nabla v\;\mathrm{d}x - \int_\Omega fv\;\mathrm{d}x - \int_{\Gamma_\mathrm{N}} gv\;\mathrm{d}S \quad + \int_{\Gamma_\mathrm{N}} \lambda v\;\mathrm{d}S + \int_{\Gamma_\mathrm{N}} \mu (u-u_\mathrm{D})\;\mathrm{d}S \quad \forall (v,\mu)\in V\times L^2(\Gamma_\mathrm{N}),$$ so we see that this modified problem admits Poisson problem $-\Delta u=f$ but with modified Neumann condition $\frac{\partial u}{\partial\mathbf{n}} = g-\lambda$ on $\Gamma_\mathrm{N}$. Simply said - requiring additional Dirichlet condition on $\Gamma_\mathrm{N}$ changed original Neumann condition required there.

If you solve this and $\lambda\ll|g|$ holds then maybe your over-specified conditions are consistent with Poisson problem. I say 'maybe' because of stability issue @ChristianClason mentioned in his answer.

Conversely we could try to enforce additional Neumann condition as a constraint (via Lagrange multiplier) on Dirichlet boundary $\Gamma_\mathrm{D}$. But that would not work either as test function $v\in V$ is zero on $\Gamma_\mathrm{D}$ which rules out constraint.

Conclusion is you that you can't expect that second order PDE will admit two independent boundary conditions.

As this question relates to FEniCS, note that it easy in FEniCS to specify $F(u,\lambda)$ and let UFL calculate its Gateaux derivative $\mathrm{D}F(u,\lambda)[v,\mu]$ by derivative() function.

Unfortunately $F(u,\lambda)$ can't be easily defined in FEniCS yet, as library does not support function $\lambda\in L^2(\Gamma_\mathrm{N})$ restricted to boundary. But another minimizations problems with Lagrange multipliers $\lambda\in L^2(\Omega)$ or $\lambda\in\mathbb{R}$ can be easily implemented this way.

• I would write minimize $\nabla ^2 u - f$ over $u\in H^1$ in the properly chosen dual space... – Jan May 14 '13 at 17:52
• My math isn't up to scratch unfortunately and I'm unsure about the mathematical implications of Banach spaces, but I struggle to see why the equation isn't a solution to a Poisson equation when the Lagrange multiplier term vanishes. From a physical point of view a solution (to the practical problem I described, I don't mean solution in the mathematical sense) must exist as far as I can see – Markus May 15 '13 at 10:52
• So it is rather an inverse problem, finding the boundary condition for the far field which, together with the Dirichlet condition that you can impose, yields the observed normal gradient at the boundary at which you measure? – Markus May 15 '13 at 12:46

Your approach cannot work, definitely because of the implementation and probably because of your formulation.

Imposing Dirichlet conditions in dolfin, eventually sets the corresponding DOFs of your testspace to zero.

This is an excerpt from the fenics-manual:

Chapter 3.3.9 (end): The application of a Dirichlet boundary condition to a linear system will identify all degrees of freedom that should be set to the given value and modify the linear system such that its solution respects the boundary condition. This is accomplished by zeroing and inserting 1 on the diagonal of the rows of the matrix corresponding to Dirichlet values, and inserting the Dirichlet value in the corresponding entry of the right-hand side vector [...]

Thus, your $v$ will be set to zero on $\Gamma_m$, what makes also your term with the multiplier disappear.

In summary, using the default routine in dolfin you cannot apply both Dirichlet and other conditions on the same boundary.

However, before you try to fix this in your implementation, go find the right test spaces for your mathematical formulation (as @Jan Blechta just has mentioned.)

• I see your point - I think my formulation might not be exactly reflecting what I've implemented - my apologies. The variational principle is but a hazy memory and I'm trying to get my head around it again. I have read the manual forward and backward together with some examples of FEniCS code implementing Lagrange multipliers. I thought the issue you raise is exactly the reason why you would use a second test function 'd'. – Markus May 15 '13 at 10:04
• I agree with @JanBlechta. At first you need to find the right space for the multiplier, which is nontrivial. Maybe texts on PDE constraint optimization, where one uses multipliers to incorporate side conditions, will give some helpful ideas. In this paper, a multiplier is used to account for time dependent Dirichlet conditions. – Jan May 15 '13 at 14:37