Simple Poisson equation in 1D with von Neumann condition on the left

Hello everyone,

I’m trying to solve the simple poisson like-equation

-u'' + u =const.

on the 1D-intervall [0,1]. Now it’s rather simple to derive the weak formulation here by just following the first couple of tutorials from the documentation. For Dirichlet-data on the left, i.e. u(0) = a and von Neumann-data on the right, u'(1) = b, it works perfectly and of course agrees with the analytical solution.
Now the problem is, that when switching the boundary conditions, i.e. u'(0) = b and u(1) =a, the numerical solution doesn’t even satisfy the von Neumann-condition on the left, i.e. at x=0. It looks, like it’s rather the negative value of this condition there, that is if I prescribe a positive slope at x=0 , the numerical solution starts with a negative slope and vice versa.I suspect, this has to do with the kind of degenerate case of surface normals to a 1D-interval.
So I was wondering what’s the recommended, proper way to deal with these scenarios?
For completeness, I’ll attach a fenics code example here:

# ------------------------------------------------------------------------
# Solve a simple non-homogeneous ODE: -d/dx(du/dx) + u = const on [0,1]
# BC: u(1) = u_D        -> Dirichlet
# BC: u'(0) = u_vN      -> von Neumann
# ------------------------------------------------------------------------

# import some useful stuff
from fenics import *
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

# create mesh and define function space
mesh = IntervalMesh(100, 0, 1)
V    = FunctionSpace(mesh, 'Lagrange', 1)

# define boundary data and righthand-side
u_D   = 1.0
u_vN  = 1.0

const = 5.0

# set up Dirichlet boundary condition
tol = 1E-14

def boundary_D(x, on_boundary):
    if on_boundary:
        if near(x, 1, tol):
            return True
            return False
        return False

bc = DirichletBC(V, u_D, boundary_D)

# define righthandside of ODE
f = Constant(const)

# now, define the variational problem:
# some infrastructure
u = TrialFunction(V)
v = TestFunction(V)

# the bilinear part
a = dot(grad(u), grad(v))*dx + u*v*dx

# the linear part
L = u_vN*v*ds + f*v*dx

# ready to solve the problem now
# define solution function
u = Function(V)

# solve the system
solve(a == L, u, bc)

# plot the solution
fsize = 15
mpl.rcParams.update({'font.size': fsize})
plt.figure(figsize=(10, 7))

x = np.linspace(0, 1, 1000)

plot(u, title="Solution", label='numerical', linewidth=3, marker='', markersize=12)
plt.legend(loc=2, fontsize=fsize)

plt.grid(b=True, which='major', color='black', linestyle=':', alpha=0.5)

Sorry, if this question has been asked before;I searched for similar topics but couldn’t find an answer.
Thanks for your help

Thats because the condition you apply is du/dn=u_N where n=-1 (from integration by parts).
Thus you need to change the sign of the boundary term.

Thanks for the fast response and I see what you mean. But isn’t this in big contrast to the same problem in 2D?
If I want to solve let’s say a Poisson-like equation in 2D and I want to modify the parts of the boundary where to apply Dirichlet- and where to apply von Neumann-data, then in the most simple case, one only has to update the ‘def boundary_D(x, on_boundary)…’-part of the code. The bilinear form itself as well as the surface part remains unchanged.
In 1D, I not only have to change the ‘def boundary’-section but additionally also the “surface-part”,which of course consists of only two single points in this case.
I thought, fenics would handle this issue somehow automatically, since the way you write down the linear and bilinear forms are very reminiscent of the multi-dimensional formulation of the PDEs of course and there nothing changes…
Anyway… thanks again for this clarification!

The difference with 1D and 2D is how you defined your boundary condition.
In 2D, its natural to set a Neumann condition as inner(grad(u),n)=du/dn=u_N, while you in your example havent included the outwards pointing normal of the domain (either -1 or +1). So in your example, you replaced du/dn*v= -u'(x)*v with u_vN*v, giving you the opposite sign.

I see; this makes it clear to me where I made my mistake…
Thanks for the reply!