DOC: uniqueness in ex05

[gdmcbain]

Not urgent. (And possibly placed somewhere better than in this issue tracker? Back on the Computational Science Stack Exchange?)

ex05 purports to solve Δ u = 0 in the unit square subject to homogeneous

• Dirichlet u = 0 on y = 0 and
• Neumann n ⋅ ∇ u = 0 on x = 0 and y = 1

conditions on three sides but driven by an integrally constrained flux on the fourth: ∫ n ⋅ ∇ u d y = 1, with the integral being over x = 1.

However, this problem, posed like this, does not have a unique solution; u [ n] = sin {(n + ½)π y} cosh {(n + ½)π x} / sinh {(n + ½)π} for n ∈ ℕ is a solution, as is any linear combination with weights summing to unity. Whereas this kind of nonuniqueness often results in a singular finite element stiffness matrix (e.g. the pure Neumann problem), this is not the case here. At first glance, the numerical solution appears to resemble u [0].

This kind of integral constraint is discussed at:

• Solving Poisson equation with current BC using FEM
• How to apply an integrated constrain condition in FEM?

The discussions imply that this is a reasonably standard boundary condition (‘I know commercial software solves this problem just as efficiently as with regular BC’) and don't mention uniqueness.

Questions:

• What is the mathematical status of this boundary condition?
• Why doesn't the finite element discretization reflect the nonuniqueness?
• Why does the discrete problem have the particular unique solution that it does?

Hunch: the discrete formulation is minimizing something which eliminates all n > 0.

[gdmcbain]

At first glance, the numerical solution appears to resemble u [0].

Not quite! It is qualitatively similar but the top-right value (x, y) = (1, 1) of the numerical solution is 0.9357… For the n-th exact solution, it is sin (n + ½)π / tanh (n + ½)π; e.g. for n = 0, coth π/2 = 1.09033…

Some plots:

Figure:— The numerical solution, ex05.x.

Figure:— The first three exact solutions, n = 0, 1, 2.

[gdmcbain]

I notice that ex05.K isn't symmetric, because ex05.B isn't. Should the integral-flux condition break symmetry?

[gdmcbain]

A bit stumped after the above, I began to wonder what this kind of boundary condition might mean. One of the posters on scicomp referred to it as an ‘electrode’ through which a given current passes. Could the implication be that the electrode is at an unknown uniform potential? Does that render the problem well posed?

[gatling-nrl]

I'm very interested in this. Did you have any more thoughts on whether the uniform potential constraint was adequate?

[gdmcbain]

'Uniform potential constraint '? Sure, but that's not in ex05. Yes, I can imagine modelling an electrode along a boundary with two constraints:

1. electrode is isopotential
2. flux through electrode matches specified value

That seems well-posed and would be easy enough to do in scikit-fem. See multipoint-constraint #438 for the isopotential (not in ex05) and use a scalar Lagrange multiplier for the flux as already in ex05.

scikit-fem/docs/examples/ex05.py

Line 54 in 71eb03e

K = scipy.sparse.bmat([[A+B, b.T], [b, None]], 'csr')

I remain puzzled by what physical problem the current ex05 solves.

[kinnala]

Currently I don't have enough time to work on this. But let's leave the issue open so maybe someone with more knowledge on the subject pops up.

[gdmcbain]

O. K.

I think I see three subproblems:

• The original problem of the uniqueness and what problem the current implementation is actually solving. This one I find to be quite a puzzle.
• What the original poster at scicomp was after and what the ‘commercial software’ does. I might ask some colleagues about that; when I first started looking at problems like ex13 in the industrial context, they were often solved with commerical software specifying a ‘current’ at a port. That might be worth an example, marginal though, in my opinion.
• Other kinds of integral constraints. The classic one is the Neumann problem, rendered unique by specifying that the solution have zero mean (e.g. Braess 2001 Finite Elements, 2nd ed., Problem III.4.21). That's probably worth an example and I expect would be implemented similarly to ex05. An interesting iterative alternative is outlined by Elman, Silvester, & Wathen 2014 Finite Elements and Fast Iterative Solvers, 2nd ed., §2.3 ‘Singular systems are not a problem’).

[kinnala]

I thought about this for a while. I wonder if this particular strong formulation is an incorrect way of looking at this. This is after all a saddle point problem so I'd think the correct strong formulation would have two unknowns, u and some Lagrange multiplier. Also, I think not all u[n] are equivalent in that they have different energies (H^1 seminorm of the solution). u[0] corresponds to the smallest energy.