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

variational formulation
1

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

2

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.

3

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!

4

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.

5

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