Catenary Problem

1

I am trying to solve the 1D Catenary Variational Problem in Fenics:

\min_y \int_{-1}^1 y \sqrt{1 + (\epsilon y')^2} \, dx

Subject to

\int_{-1}^1 \sqrt{1 + (\epsilon y')^2} \, dx = L

I am following the example here https://fenicsproject.org/qa/9024/minimization/ and my code is below. I am unable to get the Newton Solver to converge even with 10000 iterations. Am I doing something wrong with the variational form?

from fenics import *
from mshr import *

#build the problem

#parameters
ϵ = 1
L = 2

#mesh
a = -1
b = 1
Nx = 10
mesh = IntervalMesh(Nx, a, b)

#function spaces
V = FunctionSpace(mesh, "DG", 0)

#build mixed function space
P1 = V.ufl_element()
Th = P1 * P1
W = FunctionSpace(mesh, Th)

#build functions
Y = Function(W)
Y.vector()[:] = 2.7
y = TrialFunction(W)
k, λ = TestFunctions(W)

u, v = split(Y)

#minimization problem
f = u * pow(1 + ϵ ** 2 * inner(grad(u), grad(v)), 0.5) * k * dx
#constraint with lagrange multiplier
g = (pow(1 + ϵ ** 2 * inner(grad(u), grad(v)), 0.5) - L) * λ * dx

#full form
F = f + g

#derivative
J = derivative(F, Y, y)

#build the variational problem
problem = NonlinearVariationalProblem(F, Y, J=J)

#solver
solver = NonlinearVariationalSolver(problem)

solver.solve()

Thanks!

2

Keep in mind that the Lagrange multiplier in the arc length constraint is a single scalar to enforce one integral constraint over the whole domain, so it should be in a FunctionSpace of type "R", not a DG0 space where it can have a different value on each element. The following code snippet seems to solve the problem correctly:

from dolfin import *
N = 100

# Length of interval
Lhat = 2 

# Length of curve; needs to be larger than Lhat to make sense
L = Constant(10)

mesh = IntervalMesh(N,-0.5*Lhat,0.5*Lhat)
cell = mesh.ufl_cell()
Ye = FiniteElement("CG",cell,1)
# The "R" type space refers to a space with 1 DoF, for a global
# constant value.  
Le = FiniteElement("R",cell,0)
V = FunctionSpace(mesh,MixedElement([Ye,Le]))

# Set up the Lagrangian for the problem:
yl = Function(V)
y,lam = split(yl)
ds_dx = sqrt(1+(y.dx(0))**2)
constraint = (ds_dx-L/Constant(Lhat))
Lagrangian = (y*ds_dx - lam*constraint)*dx
dLagrangian = derivative(Lagrangian,yl)
bc = DirichletBC(V.sub(0),Constant(0),"on_boundary")

# Solver needs some kind of guess to set the sign of the solution.
# If you change the sign of the guess, it converges to a cantenary facing
# the other direction.  May need to adjust amplitude for different values
# of L.
initGuessAmp = 1.0
yl.interpolate(Expression(("amp*(x[0]+0.5*Lhat)*(x[0]-0.5*Lhat)","0"),
                          degree=2,Lhat=Lhat,amp=initGuessAmp))

# Need a more robust nonlinear solver than Newton iteration; try
# PETSc SNES with back-tracking ("bt") line search:
problem = NonlinearVariationalProblem(dLagrangian,yl,bc,
                                      derivative(dLagrangian,yl))
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['nonlinear_solver'] = 'snes'
prm['snes_solver']['line_search'] = 'bt'
prm['snes_solver']['maximum_iterations'] = 1000
solver.solve()

# Sanity check; should be zero
print(assemble(ds_dx*dx) - float(L))

# Plot result:
from matplotlib import pyplot as plt
plot(y)
plt.show()
1 Like