Skip to content

How to solve calculus of variations problems numerically?

An answer to this question on the Scientific Computing Stack Exchange.

Question

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?

Answer

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.