I'm trying to solve a problem of finding the displacements of an elastic material subjected to external forces. Those external forces are themselves a nonlinear function of the material displacements.

The first step is to define a functional for the energy in terms of the displacements. The solution is then the displacement which minimizes that functional over sufficiently smooth functions.

From this functional, Euler-Lagrange equations can be derived, which must be satisfied by the solution. Those equations are a system of elliptic PDEs.

To solve this problem in practice, it is necessary to discretize. This can be done directly for the functional, so that the integral of the Lagrangian over the spatial domain becomes a sum over a spatial grid, with the appropriate translation of spatial derivatives to finite differences. The problem is to minimize the sum, which is generally non-convex. The alternative is to discretize the Euler-Lagrange equations, which results in a large nonlinear system of equations. The functional contains terms which are fourth order polynomials of the displacements, while the Euler-Lagrange equations have terms which are third order polynomials of the displacements. There are no discontinuities in the problem.

My question is about the relative advantages and disadvantages of the two approaches. I have been trying quite hard to solve the second problem, i.e. solve the large nonlinear system, and I just can't seem to make progress. I have an analytic expression for the Jacobian, and I can ensure that I have a good initial guess in the sense that it is "close" to the solution as measured by the norm of the difference, but I have no guarantee that, for example, Newton's method will converge from this initial guess. I've attempted some simple globalization strategies, such as a parametrized continuation, pseudo-transient continuation, and some ad hoc regularization. The problem with Newton's method without a continuation is that it seems to always find a point at which the Jacobian becomes singular. The parametrized continuation seems to avoid this, but the Newton iterations stagnate. Quasi-Newton methods like L-BFGS work very well up to a point, and then fail catastrophically, with massive increase in the residual norm in a single iteration.

I have not attempted recasting the problem as a minimization for two reasons - first, I don't understand what advantages that would have over solving the nonlinear PDEs, and second, I would have to do a lot of calculations arising from discretizing at an early step. I realize this question may be a bit vague but I'd really appreciate any guidance - thanks!

Edit: I'll try to include some more information. One advantage I can see of discretizing the energy functional and approaching the problem as a minimization is that the Hessian would be symmetric. This contrasts with the discretization of the Euler-Lagrange equations, for which the Jacobian is not symmetric. In both approaches the Hessian/Jacobian is not positive definite. It should be possible to calculate a value which can be added to the diagonal to make a positive definite matrix, as discussed in Nocedal and Wright, but I'm unsure how to actually implement that. I suppose the idea is to put in a positive definite approximation to the Jacobian without altering the function itself.

Edit 2: Another thought about why I seem to be having so much trouble is that because my external data comes on a Cartesian grid, I am using finite difference instead of finite element. I realize that the vast majority of work in elasticity uses finite element, but I have yet to see an explicit discussion of the advantages in terms of convergence. I use centered difference operators on a 27 point stencil. The boundary conditions are homogeneous Dirichlet.