Lande-Emden type of equation and weak form

variational formulation
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