FS = 'fontsize'; LW = 'linewidth';

Suppose we wish to solve a classical problem from the calculus of variations of the form

$$ \min_{y(t)} J[y] \equiv \int_{-1}^1 f(t,y,y') dt $$

subject to the boundary conditions

$$ y(-1) = \alpha, \quad y(1) = \beta. $$

Problems like this arise in physics, for example, when $J$ is an energy functional.

A standard example of a variational problem is the catenary problem, which is to determine the shape of a hanging rope. In this case $J$ is an energy functional taking the form [1]

$$ J[y] = \int f(t,y,y') dt = \int y \sqrt{1 + (y')^2} dt. $$

To solve this problem with Chebfun, we begin with the functional and boundary conditions. We will choose boundary conditions that make it easy to verify the known solution.

f = @(y,yp) y .* sqrt(1 + yp.^2);
J = @(y) sum( f(y, diff(y)) );

alpha = cosh(-1);
beta  = cosh(1);

Since the problem is nonlinear, we will need to solve it iteratively. As an initial guess to the iteration we'll choose a simple line, and note its functional value.

dom = [-1,1];
y0 = chebfun([alpha, beta]', dom);
startJ = J(y0)
startJ =
   3.086161269630488

The simplest way to compute extrema of the functional $J[y]$ is to find solutions of the Euler-Lagrange equation

$$ \frac{\partial f}{\partial u} - \frac{d}{dx}\frac{\partial f}{\partial (u')} = 0. $$

To do so we must compute the first variations of the function $f$,

$$ \frac{\partial f}{\partial u}, \quad \frac{\partial f}{\partial u'}. $$

Using only information about first variations, we could implement a steepest descent method. However, if can also get our hands on the second variations of $f$, then we could implement a more efficient Newton method. The second- variation equivalent of the Euler-Lagrange equation, sometimes called Jacobi's accessory equation, is

$$ \frac{d}{dx}\Big(\frac{\partial^2 f}{\partial(u')^2}\Big) + \Big(\frac{d}{dx}\frac{\partial^2 f}{\partial u \partial u'}- \frac{\partial^2 f}{\partial u^2}\Big)u = 0. $$

To find the first and second variations, we will use Chebfun2. First we make a chebfun2 of $f$.

dom2 = [1 2 -2 4];
F = chebfun2(f, dom2);

Now we can compute the first and second variations of $f(y,y')$ by taking derivatives of F: each diff(F,[a b]) is a chebfun2 corresponding to

$$ \frac{\partial^{(a+b)} f}{\partial^a (u) \partial^b (u')}. $$

F1  = diff(F, [1 0]);
F2  = diff(F, [0 1]);
F11 = diff(F, [2 0]);
F12 = diff(F, [1 1]);
F22 = diff(F, [0 2]);

Now we are able to compute, say $\partial f(y,y')/\partial u$ for a particular function $y$ by writing F1(y,diff(y)) for a chebfun y.

The iteration will take Newton steps by solving the accessory equation. Five iterations are enough to get an accurate solution.

y = y0;
pref = cheboppref;
pref.plotting = 'off';
pref.display = 'none';

for k = 1:5
    % The first and second variations of f.
    yp = diff(y);
    f1  = F1(y, yp);
    f2  = F2(y, yp);
    f11 = F11(y, yp);
    f12 = F12(y, yp);
    f22 = F22(y, yp);

    % Compute the next Newton step by solving the accessory equation.
    N = chebop( @(x,u) diff(f22.*diff(u)) + (diff(f12)-f11).*u, dom );
    N.lbc = 0;
    N.rbc = 0;
    u = solvebvp(N, (f1-diff(f2)), pref);

    % Increment by the step.
    y = y + u;

    % Print out the current value of the energy functional.
    nextJ = J(y)
end
nextJ =
   2.840828750691852
nextJ =
   2.816226754422263
nextJ =
   2.813551498179433
nextJ =
   2.813430779819845
nextJ =
   2.813430203941191

The exact solution is a simple $\cosh$. Here is the error in our solution.

y_exact = chebfun('cosh(x)', dom);
norm(y - y_exact)
ans =
     7.515155564171490e-06

And the values of the energy functional:

fprintf('\n  final J[y]: %.16f\n', J(y))
fprintf('optimal J[y]: %.16f\n', J(y_exact))
  final J[y]: 2.8134302039411909
optimal J[y]: 2.8134302039235082

Finally, a plot of the catenary:

plot(y_exact, 'r--', LW,2), hold on
plot(y, 'k-', LW,1)
axis equal
title('Solution to the catenary problem', FS, 12)
xlabel('x', FS, 14), ylabel('y(x)', FS, 14)
legend('exact','computed','location','north')

References

  1. Charles Fox, An Introduction to the Calculus of Variations, Courier Dover Publications, 1987.