Application Center - Maplesoft

App Preview:

The VariationalCalculus Package

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


The Calculus of Variations in Maple 8

Copyright 2002 Waterloo Maple Inc


The new VariationalCalculus  package provides routines for solving problems in the calculus of variations, which studies nature's most "efficient" curves and surfaces.  Classical problems from the calculus of variations include:

  • Find the shortest path between two points on a 3-D surface, such as a cone (a geodesic problem).
  • Shape a ramp between two heights such that a ball rolling down it reaches the bottom in minimum time (the Brachistochrone problem).
  • Find the shape of a soap film having minimum surface area spanning a given wire frame.

Such problems can often be solved with the Euler-Lagrange equation, which generalizes the Lagrange Multiplier Theorem for minimizing functions of real variables subject to constraints.  The Euler-Lagrange equation is easy to write down in general but notoriously difficult to write down and solve for most practical problems.  The VariationalCalculus  package automates the construction and analysis of the Euler-Lagrange equation.


Example: The Brachistochrone Problem

The Brachistochrone problem can be stated as follows: Given two endpoints in the plane, find the curve y(x) between them such that a ball of unit mass rolls along the curve under the influence of gravity in minimum time.

>    restart;

>    with(VariationalCalculus);

[ConjugateEquation, Convex, EulerLagrange, Jacobi, Weierstrass]

First we write down the falling time over an infintesimal distance dx in terms of y(x)  and yInit , assuming the gravitational constant is 1.  This is found in standard textbooks on classical mechanics.

>    fallTime := sqrt( (1+diff(y(x),x)^2)/(2*(yInit-y(x))) );

fallTime := ((1+diff(y(x),x)^2)/(2*yInit-2*y(x)))^(1/2)

We then use the EulerLagrange  function to compute the Euler-Lagrange equations for this functional in terms of y(x)  and its derivatives.  This returns a set of ODEs.

>    EL := EulerLagrange( fallTime, x, y(x) ) ;

EL := {((1+diff(y(x),x)^2)/(2*yInit-2*y(x)))^(1/2)-diff(y(x),x)^2/((1+diff(y(x),x)^2)/(2*yInit-2*y(x)))^(1/2)/(2*yInit-2*y(x)) = K[1], 1/((1+diff(y(x),x)^2)/(2*yInit-2*y(x)))^(1/2)*(1+diff(y(x),x)^2)/(...
EL := {((1+diff(y(x),x)^2)/(2*yInit-2*y(x)))^(1/2)-diff(y(x),x)^2/((1+diff(y(x),x)^2)/(2*yInit-2*y(x)))^(1/2)/(2*yInit-2*y(x)) = K[1], 1/((1+diff(y(x),x)^2)/(2*yInit-2*y(x)))^(1/2)*(1+diff(y(x),x)^2)/(...
EL := {((1+diff(y(x),x)^2)/(2*yInit-2*y(x)))^(1/2)-diff(y(x),x)^2/((1+diff(y(x),x)^2)/(2*yInit-2*y(x)))^(1/2)/(2*yInit-2*y(x)) = K[1], 1/((1+diff(y(x),x)^2)/(2*yInit-2*y(x)))^(1/2)*(1+diff(y(x),x)^2)/(...

Extract the first ODE from the list

>    ODE1 := remove(has, EL, diff(y(x),x,x))[1];

ODE1 := ((1+diff(y(x),x)^2)/(2*yInit-2*y(x)))^(1/2)-diff(y(x),x)^2/((1+diff(y(x),x)^2)/(2*yInit-2*y(x)))^(1/2)/(2*yInit-2*y(x)) = K[1]

Compute y(x)  as a parametric curve (x(s), y(s)) using the dsolve  command with the parametric  option.

>    sol := dsolve(ODE1, parametric, y(x));

sol := [y(_T) = 1/2*(2*yInit*K[1]^2*_T^2-1+2*yInit*K[1]^2)/(1+_T^2)/K[1]^2, x(_T) = 1/2*(_T+arctan(_T)+arctan(_T)*_T^2+2*_C1*K[1]^2+2*_C1*K[1]^2*_T^2)/(1+_T^2)/K[1]^2]

>    sol := subs( {_T=s}, sol);

sol := [y(s) = 1/2*(2*yInit*K[1]^2*s^2-1+2*yInit*K[1]^2)/(1+s^2)/K[1]^2, x(s) = 1/2*(s+arctan(s)+arctan(s)*s^2+2*_C1*K[1]^2+2*_C1*K[1]^2*s^2)/(1+s^2)/K[1]^2]

>    assign(sol);

Specify the endpoints of the curve: (0,1)  and (1, 1/2)

>    xInit := 0; yInit := 1; xEnd := 1; yEnd := 1/2;

xInit := 0

yInit := 1

xEnd := 1

yEnd := 1/2

The parameter s  is not equivalent to time t , so we must solve for the value of s  at which y(s) = yEnd

>    solve(y(s)=yEnd, s);

(1-K[1]^2)^(1/2)/K[1], -(1-K[1]^2)^(1/2)/K[1]

>    sEnd := %[1];

sEnd := (1-K[1]^2)^(1/2)/K[1]

Use the endpoints to solve for the unknown constants K[1]  and _C1

>    _C1 := solve(subs({s=sEnd}, x(s)) = xEnd, _C1);

_C1 := 1/2*(-(1-K[1]^2)^(1/2)*K[1]-arctan((1-K[1]^2)^(1/2)/K[1])+2*K[1]^2)/K[1]^2

The finite-length curve (x(s), y(s))  is parametrized over [-infinity, infinity] , so xInit = x(-infinity)

>    limit(x(s), s=-infinity);


>    K[1] := fsolve( % = xInit, K[1]);

K[1] := .9832314847

The path along which a ball would roll in minimum time under the influence of gravity.  The result is surprising because the curve actually dips below the right end point.

>    plot([x(s),y(s), s=-10..10], xInit..xEnd, 0..yInit);

[Maple Plot]