# 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? Oct 22, 2019 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. Oct 22, 2019 at 23:47
• You may want to take a look at this paper [PDF]. Oct 23, 2019 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. Oct 23, 2019 at 8:27

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] == 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?