# How to solve calculus of variations problems numerically?

For example, how to solve the well-known isoperimetric problem (i.e., to enclose the largest area with a fixed-length curve)?

We can simplify things a bit and fix the two ends of the curve at $$[a,0]$$, $$[b,0]$$, then the problems is

$$\text{maximize} \quad \int_a^b y(x) dx \quad \text{subject to} \quad \int_a^b \sqrt{1+y’(x)^2} dx = L$$

How to directly solve it without applying the Euler-Lagrange equation? It's a classic optimization problem, right? A quick search leads to Euler's finite difference method and Ritz's method, and some examples with fixed-end constraint. So I guess the fixed-length constraint should be handled via Lagrangian multiplier.

Are there "better", more state-of-the-art method? Also are there numerical packages that can take an integral objective and some arbitrary constraints, and give you the solution without much user intervention?

• "How to directly solve it without applying the Euler-Lagrange equation?" I think that's somewhat impossible... What's wrong with Euler-Lagrange approach? – Alone Programmer Oct 22 '19 at 22:08
• @AloneProgrammer how dare I to say something is wrong with Euler/Lagrange's work :) I am asking because first, the resultant EL equation could be too complex to solve--you have to rely on numerics anyway; second, to me we already have a well-defined optimization problem, we shouldn't have to simplify any further just to solve it. – Taozi Oct 22 '19 at 23:47
• You may want to take a look at this paper [PDF]. – Rodrigo de Azevedo Oct 23 '19 at 7:25
• I believe we can make the $=$ in the constraint a $\leq$ without changing the optimal value. If we do that, this problem becomes convex (and remains so when we discretize into finite dimensions). In particular I believe this turns into a second order cone program. Can look into the details more tomorrow maybe. – cdipaolo Oct 23 '19 at 8:27

## 2 Answers

In the specific problem you ask about (unlike Richard's more general answer), it turns out that you can relax the $$=$$ constraint into a $$\leq$$ without changing the optimal value, and that the resulting convex problem can be solved with the CVX software Richard mentioned.

Details: Intuitively the relaxation is possible since if the function had arc length strictly less than $$L$$, you could scale it up to have arc length $$L$$ exactly, but scaling the function up will necessarily increase the integral. As a result, if we solve the optimization problem with a $$\leq$$ replacing the $$=$$ in the constraint, the optimal value and function $$y(x)$$ is the same.

It turns out the resulting problem is convex in the parameter $$y(x)$$, albeit infinite dimensional. The objective is linear. Moreover, the function $$\sqrt{1 + y'(z)^2} = \|(1, y'(z))\|_2$$ is a convex function of $$y(x)$$ as a norm of a linear function of $$y(x)$$ (the derivative is linear). Thus, since the integral of a family of convex functionals is convex, the mapping $$y(x) \mapsto \int_a^b \sqrt{1 + y'(x)^2}\,dx$$ is convex. The constraint is then convex, as it asks for $$y(x)$$ to lie in a sub-level set of a convex function of $$y(x)$$.

The following Python code gives an example of solving this problem with the CVXPY software Richard mentioned:

import numpy as np
import cvxpy as cp

# Problem Data
n = 101 # points in interval [0,1] to discretize y at
L = 1.5

y = cp.Variable(n) # y[i] = y(i/n) for i=0,1,...,n
x = np.arange(n) / n
dy = y[1:] - y[:-1]
dx = x[1:] - x[:-1]

# Construct Optimization Problem
objective = (cp.sum(y[:-1]) + cp.sum(y[1:])) / 2 / n
constraints = [
y == 0, y[-1] == 0,
cp.sum(cp.norm2(cp.vstack([dy, dx]), axis=0)) <= L,
]

prob = cp.Problem(cp.Maximize(objective), constraints)

# Solve Problem
prob.solve(verbose=True)


The solution (found at y.value) is given as follows: As Richard points out, you can't hope to be able to add generic constraints to this problem and continue to have access to a simple and/or efficient solver. Nevertheless, you can add convex constraints like this arc-length one in a pretty straightforward way without too much trouble. For example, the following gives the output of the same optimization problem if I now add the constraint $$y(1/2) \geq 0.55$$: I think the meat of your question is:

Are there numerical packages that can take an integral objective and some arbitrary constraints and give you a solution?

The answer is: mostly, no.

However, a few packages, such as GPOPS, JuMP, and PyOMO, have facilities for this which sometimes work if your problem is of the right form.

For particular classes of problems there are excellent solvers, but arbitrary constraints allow you to actually build quite a bit of complexity into a problem pretty easily. Consider, for instance singular control problems.

For instance, if you constrain a variable to belong to a set, then your problem becomes an integer program. Those are in the NP complexity class, so we don't and probably won't have good, general-purpose algorithms that solve all instances.

If the function you're integrating curves the wrong way you move from having a convex program (usually solvable easily) to a non-convex program, where, again, there's not a general algorithm. Or, if you have a case-based function, then you may no longer have smooth or continuous derivatives.

Even if your problem happens to fall into a class of instances for which we have good algorithms, you have to explain that to the computer. Packages like cvxpy require you to write your problem using a disciplined library of functions which not only restrict you to solvable instances, but which actually constitute a recipe for the solver to follow.