How to solve 1d equation with Integral and Robin bcs

variational formulation
1

Hello everyone!

I am curious how to solve the following problem using Fenics

\dfrac{d^2 u}{d x^2} - u = 0, \quad x \in [0,1],

with the following BCs

\int_0^1 u\; dx = 1, \quad \dfrac{du}{dx}=u |_{x=1}

It’s obvious that it has a solution u_{ex} = \dfrac{1}{e-1}exp(x), but i’am struggling with variational formulation for this.

I’ve seen approach for pure Neumann conditions with Lagrange multipliers to include integral condition.
However, here it’s slightly different, due to the fact that we have Robin type BC only on one part of the boundary.

I came up with smt like this, but it gives wrong results
u v |_{x=1} - \dfrac{du}{dx} v |_{x=0} - \dfrac{du}{dx} \dfrac{dv}{dx} dx - u v \; dx + \lambda v \; dx + test(\lambda)(u-1) \; dx = 0,
where \lambda is a Lagrange multiplier.

Here is the code

from dolfin import *
import math as mt

# definig mesh
R_out = 1
R_in = 0
h = R_out-R_in
mesh = IntervalMesh(20, R_in, R_out)
tol = 1E-14

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], R_in)
    
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], R_out)    

r_ = Expression('x[0]',degree=2)

left = Left()
right = Right()
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
ds = Measure("ds", subdomain_data=boundaries)


P2 = FiniteElement("Lagrange", interval, 1)
R = FiniteElement("R", interval, 0)
element = MixedElement([P2, R])
V0 = FunctionSpace(mesh, element)

U_r = TrialFunction(V0)
V_r = TestFunction(V0)

trial_Ur, trial_lam = split(U_r)
test_Vr, test_lam = split(V_r)

# the bilinear part
a = ( -trial_Ur.dx(0)*test_Vr.dx(0) -  trial_Ur*test_Vr )*dx 
a += ( test_lam*(trial_Ur) + trial_lam*test_Vr)*dx 
a += - trial_Ur.dx(0)*test_Vr*ds(1) + trial_Ur*test_Vr*ds(2) 
# the linear part
L= (test_lam/h)*dx

u = Function(V0)
# solve the system
solve(a == L, u)
2

Ok, guys. For those of you who will be interested in solution
I’ve found the way to implement it using a little bit different approach, without Lagrange multipliers

Now, we consider the following equivalent problem
\dfrac{d^2 u}{d x^2} - (u+c) =0 , \quad x \in [0,1]
with BCs like this
\int (u+c) dx = 0, \quad \dfrac{du}{dx} = (u+c)|_{x=1}, \quad u|_{x=0}=0.

Thus, we derive the weak formulation

(u+c)v|_{x=1} - \dfrac{d u}{d x} \dfrac{d v}{d x} dx - (u+c) v dx + (u+c) \; test(c) dx =0

The solution we are interested in is obtained by summation of u and c.

Now it works

from dolfin import *
import math as mt

# definig mesh
R_out = 1
R_in = 0#0.5
h = R_out-R_in
mesh = IntervalMesh(10, R_in, R_out)
tol = 1E-14

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], R_in)
    
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], R_out)    

r_ = Expression('x[0]',degree=2)

left = Left()
right = Right()
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
ds = Measure("ds", subdomain_data=boundaries)



# definig Function space on this mesh using Lagrange polynoimals of degree 2.
P2 = FiniteElement("Lagrange", interval, 2)
R = FiniteElement("R", interval, 0)
element = MixedElement([P2, R])
V0 = FunctionSpace(mesh, element)


f = Expression ( "0.0", degree = 0 )

U_r = TrialFunction(V0)
V_r = TestFunction(V0)

trial_Ur, trial_lam = split(U_r)
test_Vr, test_lam = split(V_r)
bcs = DirichletBC(V0.sub(0), 0.0, boundaries, 1)

# the bilinear part
a = ( -trial_Ur.dx(0)*test_Vr.dx(0) -(trial_Ur+trial_lam)*test_Vr )*dx 
a += ( test_lam*(trial_Ur + trial_lam))*dx + (trial_Ur+trial_lam)*test_Vr*ds(2)  #- trial_Ur.dx(0)*test_Vr*ds(1)
# the linear part
L=f*test_Vr*dx + (test_lam/h)*dx

u = Function(V0)
# solve the system
solve(a == L, u, bcs)