share/calculus/optvar.dem

   1 load('optvar)$
   2 /* This macsyma batch file illustrates how to use the
   3 functions in OPTVAR > or OPTVAR LISP to help solve classical
   4 variational or optimal control problems.  These functions
   5 use the calculus of variations and Pontryagin's
   6 Maximum-Principle to derive the governing differential and
   7 algebraic equations; and this demonstration also shows how
   8 the governing equations may be solved using the built-in
   9 SOLVE function and/or the functions in the SHARE files  ODE2
  10 LISP, and DESOLN LISP.  For more details, see the
  11 corresponding SHARE text files OPTVAR USAGE, ODE2 USAGE, or
  12 DESOLN USAGE.
  13
  14 The illustrative examples here are intentionally elementary.
  15 For a more thorough discussion of the mathematical
  16 principles demonstrated here, and for results with more
  17 difficult examples, see the the ALOHA project technical
  18 report by David Stoutemyer, "Computer algebraic manipulation
  19 for the calculus of variations, the maximum principle, and
  20 automatic control", University of Hawaii, 1974.
  21
  22 Many classical variational problems are analogous to special
  23 cases of choosing a path which minimizes the transit time
  24 through a region for which the speed is a function of
  25 position.  There are applications in, optics, acoustics,
  26 hydrodynamics, and routing of aircraft or ships. In the
  27 two-dimensional case, assuming the path may be represented
  28 by a single-valued function, Y(X), the transit time between
  29 Y(A) and Y(B) is given by the integral of
  30 Q(X,Y)*SQRT(1+'D(y,x)**2), from X=A to X=B, where D is an
  31 alias for DIFF and Q(X,Y) is the reciprocal of the speed,
  32 and where we have used the fact that 'D(arclength,X) =
  33 SQRT(1+'D(Y,X)**2).  For simplicity, assume Q independent of
  34 X.  We may use the function EL, previously loaded by
  35 LOADFILE(OPTVAR,LISP,DSK, SHARE), to derive the associated
  36 Euler-Lagrange equation together with an associated energy
  37 and/or momentum integral if they exist: */
  38 aa:el(q(y)*sqrt(1+'d(y,x)**2), y, x);
  39 /* To expand the Euler-Lagrange equation: */
  40 dependencies(y(x))$
  41 part(aa,2), eval, d;
  42 assume(k[0]>0)$
  43 /* The routines in the SHARE file ODE2 LISP, previously
  44 loaded, may solve an expanded first or second order
  45 quasi-linear ordinary differential equations; so the
  46 equation must be linear in its highest- order derivative.
  47 The Euler-Lagrange equation is always of this form, but when
  48 given a second-order equation, the ODE2 solver often returns
  49 with a first-order equation which we must quasi-linearize
  50 before proceeding; so it is usually most efficient to take
  51 advantage of a first integral when one exists, even though
  52 it requires a certain amount of manipulation.  The SOLVE
  53 function is currently somewhat weak with fractional powers;
  54 so we must massage the above energy integral before solving
  55 for 'D(Y,X): */
  56 part(aa,1) * sqrt(1+'d(y,x)**2), expand, eval;
  57 solve(%/k[0], 'd(y,x));
  58 ode2(ev(part(%,2),eval), y, x);
  59 /* We must specify Q(Y) to proceed further.  Q(Y)=1 is
  60 associat- ed with the Euclidian shortest path (a straight
  61 line), Q(Y)=SQRT(Y) is associated with the stationary
  62 Jacobian-action path of a projectile (a parabola), Q(Y)=Y is
  63 associated with the minimum-surface body of revolution (a
  64 hyperbolic cosine), and Q(Y)=1/SQRT(Y) is associated with
  65 the minimum-time path of a falling body, starting at Y=0,
  66 measured down, (a cycloid).  Of these, the cycloid presents
  67 the most difficult integration.  In fact, I have never seen
  68 the the integration for this case performed directly except
  69 by macsyma: */
  70 ratsubst(1/sqrt(y),q(y),%);
  71 %,integrate;
  72 /* This equation may be simplified by combining the
  73 LOGs and clearing some fractions, but it is transcendental
  74 in Y; so there is no hope for a closed-form representation
  75 for Y(X).  For completeness we should try the other
  76 alternative for 'D(Y,X), but it turns out to lead to the
  77 same result.
  78
  79 Most authors solve this problem by introducing the change of
  80 variable 'D(Y,X)=TAN(T), useful also for other Q(Y), which
  81 leads to the para- metric representation: Y=K[0]*COS(T)**2,
  82 X=C+K[0]*(T+SIN(T)*COS(T)). Solving the former for T and
  83 substituting into the latter gives the representation
  84 equivalent to that found directly by macsyma.
  85
  86 Some straightforward investigation reveals that  C  is the
  87 initial value of  X,  but the equation is transcendental in
  88 K[0], so K[0] cannot be determined analytically.  However,
  89 it may be determined to arbitrary numerical accuracy by an
  90 iterative algorithm such as that described in the SHARE file
  91 ZEROIN USAGE.
  92
  93 To prove that the above solution is a strong global minimum,
  94 we must prove that such a minimum exists and that no other
  95 answer gives a smaller value to the functional.  This may be
  96 done by the direct approach of Caratheodory,  or by checking
  97 the classical Weierstrass, Weierstrass-Erdman, Legendre, and
  98 and Jacobi necessary and sufficient conditions.  Computer
  99 algebraic manipulation can assist these investigations, but
 100 they are beyond the scope of this demonstration.  Instead,
 101 to keep the demonstration brief, we will now consider other
 102 kinds of problems.
 103
 104 Stationary values of a functional may be subject to
 105 algebraic or differential constraints of the form
 106 F(X,Y(X),'D(Y,X))=0 and/or isoperimetric constraints of the
 107 form INTEGRATE(F(X,Y(X),'D(Y,X)),X,A,B) = J, where J is a
 108 given constant. Sometimes it is possible to eliminate one or
 109 more constraints, substituting them into the functional and
 110 any remaining constraints; but if not, we may use Lagrange
 111 multipliers.  For example, suppose we wish to determine the
 112 curve that a string of given length assumes to minimize its
 113 potential energy.  The potential energy is
 114 INTEGRATE(Y*SQRT(1+'D(Y,X)**2), X, A, B), whereas the length
 115 is INTEGRATE(SQRT(1+'D(Y,X)**2), X, A, B); so letting MULT
 116 be the Lagrange multiplier, our augmented functional is of
 117 the form that we have been studying, with Q(Y)=Y+MULT.
 118 Consequently, we may simply substitute this in the general
 119 integral that we found before: */
 120 subst(q(y)=Y+mult, %th(3));
 121 /* Type  NONZERO;  in response to the question
 122 generated by the following statement: */
 123 %,integrate;
 124 /* to get Y as a function of X: */
 125 %/k[0];
 126 exp(lhs(%))=exp(rhs(%));
 127 solve((%-2*y-2*mult)**2, y);
 128 %,eval,expand,radexpand:true;
 129 /* We may rewrite the above expression as: */
 130 k[0]*cosh(x/k[0]+c) - mult $
 131 /* We may substitute this into the integral constraint
 132 to get 1 equation relating the constant K[0] and C to the
 133 given length L: */
 134 sqrt(1+d(%,x)**2);
 135 /* This is simply COSH(X/K[0]+C), so: */
 136 l = integrate(cosh(x/k[0]+c), x, a, b);
 137 /* Given A, B, Y(A), and Y(B), we may substitute into
 138 the general solution to get two more relations among K[0], C
 139 and MULT, but the three equations are transcendental; so a
 140 numerical solution is generally necessary.
 141
 142 For completeness, we should also investigate the case of
 143 K[0]=0, which yields the solution  Y=0; but for brevity, we
 144 will omit that analysis. */
 145 "end of part 1"$