Lande-Emden type of equation and weak form

variational formulation
1

Hello everyone,
I am hoping for some help concerning PDEs, where the boundary conditions (Dirichlet and von Neumann) are specified for the very same boundary.
The Lane-Emden equation is an ODE governing hydrostatic equilibrium for simple gases in spherical symmetry but for the sake of simplicity let’s use something like

u'' - u = 0

as an example. Now prescribing boundary conditions like u(0) = 0 and u'(1)=1, it’s easy to derive the weak form and implement it in Fenics.
However, boundary conditions like u(0)=0 and u'(0)=1 are also possible and the analytic solution to this BVP is easily written down.
This type of boundary conditions also appear in the Lane-Emden equation; i.e. the value of the solution and its derivative at the same point are prescribed.
Here, I am struggling to arrive at the proper weak formulation where one typically applies Dirichlet- and von Neumann-conditions to different parts of the boundary.
Of course, one cheap trick would be to neglect the von Neumann-condition in the first step and iteratively change an additional Dirichlet-condition at x=1 while monitoring the first derivative of the numerical solution at x=0 and adjusting the “fake” Dirichlet condition at the right boundary point accordingly.
However, this is not a very elegant solution and rather inferior to let’s say a Finite-Difference-implementation, where one simply “integrates outwards”, starting from x=0.
So I was wondering whether I’m missing some obvious step and I would appreciate any helpful comments…

Thanks very much!

2

Just thinking off the top of my head, you could try a Nitsche-like formulation, although my quick test only converges optimally in H^1, without attaining a higher rate in L^2, even with the analog to the “symmetric” Nitsche method variant (which is no longer symmetric when modified to apply two BCs to one end of the domain and none to the other):

from dolfin import *
N = 32
mesh = UnitIntervalMesh(N)
V = FunctionSpace(mesh,"CG",1)
u = TrialFunction(V)
v = TestFunction(V)

# Manufacture exact solution to check accuracy:
x = SpatialCoordinate(mesh)[0]
u_ex = (1.0/pi)*sin(pi*x)
f = -div(grad(u_ex)) + u_ex # Generalize to nonzero volume source term
g = Constant(0.0) # Dirichlet data at left endpoint
h = Constant((1.0,)) # Neumann data at left endpoint

# Characteristic functions for endpoints when restricted to ds measure
leftChar = 1.0-x
rightChar = x
n = FacetNormal(mesh)

# PDE:
interiorRes = dot(grad(u),grad(v))*dx + u*v*dx - f*v*dx

# Weak BCs:  Idea is to keep the boundary term from integration-by-parts at
# the right endpoint, while applying Nitsche's method at the left endpoint,
# but modifying the consistency term to enforce a Neumann BC.
he = CellDiameter(mesh)
pen = Constant(1e0)/he
gamma = Constant(1.0) # (Change to -1 for non-symm, which works with pen=0)
bdryRes =  (leftChar*(-dot(h,n)*v - gamma*(u-g)*dot(grad(v),n) + pen*(u-g)*v)
            + rightChar*(-dot(grad(u),n)*v))*ds

# Solve and check error:
res = interiorRes + bdryRes
uh = Function(V)
solve(lhs(res)==rhs(res),uh)

# Optimal:
print("H1 error = "+str(sqrt(assemble((grad(uh-u_ex))**2*dx))))

# Sub-optimal:
print("L2 error = "+str(sqrt(assemble((uh-u_ex)**2*dx))))
1 Like
3

Thanks for the reply, I will take a look at it!
I didn’t imagine that it would be so elaborate for a rather elementary problem, in particular when compared to the trivial implementation in e.g. finite differences…