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
else:
return False
else:
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.xlabel('x')
plt.ylabel('u(x)')
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