Section 47-02.html

Unit 9: Calculus of Variations

Chapter 47: Constrained Optimization

Section 47.2: Queen Dido's problem

Copyright

All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior written permission of the publisher. Printed in the United States of America.

Initializations

> restart;

> libname:="C:/Program Files/Maple 6.01/Calcvar",libname:

> with(Calcvar):

> libname:="C:/Program Files/Maple 6.01/Partials",libname:

> with(Partials):

> with(student):
with(plots):

Warning, the name changecoords has been redefined

>

The Isoperimetric Problem

The problem of stationarizing the functional

subject to the condition that should satisfy

where is a fixed constant, is called the isoperimetric problem because the constraint demands that have a fixed length. (The constraint integral is simply the arclength for .) The word isoperimetric is a derivation of the Greek words for same and perimeter .

A most amusing version of this problem is given in [1]. Supposedly, in 850 B.C., Queen Dido of Carthage, a refugee from Tyria, asked the North African chieftain named Jarbas for as much land as could encompassed with the hide of a cow. Queen Dido's solution was to cut the hide into strips, and to lay the strips out in a semicircle with the coast of North Africa as a diameter. By this account, the enclosed area was maximized, and the city-state of Carthage founded.

Queen Dido solved the isoperimetric problem with free endpoints by enclosing the maximum area with a semicircle. The isoperimetric problem with fixed endpoints, the version of the problem we solve first in this section, does not turn out to be a semicircle! However, we still refer to it as Queen Dido's problem with fixed endpoints.

>

Queen Dido's Problem with Fixed Endpoints

In the fixed-endpoint version of Queen Dido's problem, we are to find the curve which has a fixed and given perimeter , and which, together with the interval on the -axis, encloses a maximum area.

The functional to maximize is

and the constraint is

(Of course, there is no loss in generality by arranging the endpoints symmetrically with respect to the origin.)

To solve this problem with the formalism of the calculus of variations, write, according to Table 47.1 of Section 47.1,

as the integrand of the modified functional . This integrand does not explicitly contain the independent variable , so from Section 46.4 we immediately have the first integral

that is,

Separation of variables gives for the solution of the Euler-Lagrange equation

the equation of a circle centered at . This is an implicit form of the function we are trying to determine. It contains three parameters, , and . Two conditions we can impose are the endpoint conditions

= 0

resulting in

and

from which we get, by subtracting the second from the first, . This tells us that because . Hence, the two endpoint equations both reduce to

= +

and

The arc length constraint them becomes

Eliminating from the arc length constraint and leads to

as an implicit definition of .

Let us next implement these calculations in Maple. To use the formalism of the Calcvar package, establish the dependence of on via

> depends(y,x);

Warning, alias or macro ddy defined in terms of y

Warning, alias or macro dy defined in terms of y

Warning, alias or macro y'' defined in terms of y

Warning, alias or macro y' defined in terms of y

>

and write

> F := y+m*sqrt(1+dy^2);

>

where we have opted to use the single rather than which requires six characters to create in Maple. This is a case where the integrand is free of the independent variable , so a first integral exists. Using Maple's Euler_Lagrange command, we obtain

> q :=Euler_Lagrange(F,x,y);

$q := {1+m*`y'`^2*diff(y,`$`(x,2))/((1+`y'`^2)^(3/2)...$

>

We have obtained both the first integral and the expanded form of the Euler-Lagrange equation. Extracting the first integral, and simplifying it somewhat, we have

> Q := normal(op(select(has,q,K[1])));

>

Let us also replace the cumbersome with the simpler , obtaining

> Q1 := subs(K[1]=k,Q);

>

This first integral can be integrated by separating variables, but we prefer to obtain the solution with Maple's dsolve command. We get

> q1 := dsolve(Q1,y);

>

Selecting one branch, and replacing the cumbersome integration constant with the simpler , we have

> q2 := subs(_C1=d,q1[1]);

>

The remaining manipulations belong to the domain of algebra. The calculus of variations is behind us now, and it is wise to have access to an unaliased variable . Hence,

> free(y,x);

>

to remove the alias attached to . However, at this point, Maple has recorded

> q2;

>

for the solution of the Euler-Lagrange equation. Hence, we must also replace with via the substitution

> q4 := subs(y(x)=y,q2);

>

It helps to remove the square root by squaring both sides. Beforehand, however, let's move the to the other side of the equation so that the radical is isolated on the left. Subtracting from both sides of the equation is the precise way to move it to the other side. Hence, we have

> q5 := q4 - (k = k);

>

Now square both sides, obtaining

> q6 := map(x -> x^2, q5);

>

Complete the square on the left is as useful as it is aesthetic, and we find

> q7 := completesquare(q6,y);

>

If we bring the term in to the right by subtracting it from both sides, we obtain

> q8 := q7 - ((y-k)^2=(y-k)^2);

>

which is the equation of a circle. This is an implicit form of the function we are trying to determine. It contains three parameters, , , and . Two conditions we can impose are the endpoint conditions = 0. This results in the two equations

> eq1 := subs(x=-a,y=0,q8);
eq2 := subs(x=a, y=0,q8);

>

from which we get, by subtracting the second from the first,

> q9 := simplify(eq1-eq2);

>

This tells us that because . Hence, the two endpoint equations both reduce to

> eq3 := subs(d=0,eq1);

>

Setting this relationship aside for the moment, we turn our attention to the integral constraint. We must apply the constraint to the solution in order to obtain a second equation relating the parameters and . The constraint equation is an arclength integral, so we need to compute . We opt for implicit differentiation implemented via Maple's implicitdiff command. Since the solution to which we will apply this command contains the parameter , we take the trouble to set that parameter to zero upon differentiating.

Hence, is

> q10 := subs(d=0,implicitdiff(q8,y,x));

>

But the expression for contains itself, so we need to solve for explicitly. We get

> q11 := subs(d=0,[solve(q8,y)]);

>

where we again made sure to set to zero. Picking a branch, and replacing with it in the expression for , we get

> q12 := subs(y=q11[2],q10);

>

The integrand of the constraint integral is , which is then

> q13 := simplify(sqrt(1+q12^2),symbolic);

>

Integrating in Maple gives

> int(q13,x=-a..a);

>

Again, Maple's propensity to choose the log form of the inverse trig functions is a disaster. Rather than grapple with this most inhospitable result, we simply consult an elementary calculus text and write the antiderivative of

> q14 := m*arcsin(x/m);

>

For sure, the derivative of this antiderivative is again , as we verify by

> simplify(diff(q14,x),symbolic);

>

Evaluating this antiderivative at the endpoints, we get

> q15 := simplify(subs(x=a,q14) - subs(x=-a,q14)) = c;

>

as the value of the constraint integral. To eliminate the parameter , we solve the relationship

> eq3;

>

for , and make the obvious substitution. However, there are two solutions for , as we see by

> q16 := solve(eq3,m);

>

so we select the first, and with it, obtain the equation

> q17 := subs(m=q16[1],q15);

>

Since and are given, this is a single equation in the single unknown . It cannot be solved explicitly for , but it can be slightly simplified if the arcsine function is isolated via

> q18 := q17/2/sqrt(a^2+k^2);

Computing the sine of each side then gives

> q19:=map(sin,q18);

>

as an implicit definition of .

>

Example 47.1

As a specific example of this version of Queen Dido's problem with fixed endpoints, take and . Then, the parameter is determined by the equation

> q20 := subs(a=1,c=2*Pi/3,q19);

>

A numeric solution is

> fsolve(q20,k);
fsolve(q20,k,-2..0);

>

which looks suspiciously like + , a guess we confirm with

> sqrt(3.);

>

Surprisingly, Maple does not recognize that the equation is just and with = 2, we are looking at . Hence, really is + .

To graph , we need to compute = + = + = + 2.

This gives implicitly as

> q21 := subs(m=2,d=0,k=-sqrt(3),q8);

>

where we have chosen so the circle has its center below the -axis. That assures us the arc of the circle, namely, , graphed in Figure 47.1, lies above the -axis.

> implicitplot(q21,x=-1..1,y=0..2,scaling=constrained, xtickmarks=3, ytickmarks=3, labels=[x,`y `], labelfont=[TIMES,ITALIC,12]);

>

To find the area enclosed by this arc, we need to obtain explicitly. Hence, solve for , obtaining the two branches

> q22 := solve(q21,y);

>

The appropriate branch is the one with the plus sign in front of the radical. Hence,

> if has(q22[1],-sqrt(4-x^2)) then Y:=q22[2] else Y:=q22[1]; fi;

>

selects the appropriate branch, and

> q23 := int(Y,x=-1..1):
q23 = evalf(q23);

>

gives the area encompassed by and that portion of the -axis in the interval .

>

Queen Dido's Problem with Free Endpoint

With fixed arclength , enclose the maximum area between the curve and the -axis, if characterizes the fixed left endpoint, and characterizes the free right endpoint. The transversal is , and must be determined by the transversality condition.

Thus, maximize the functional

subject to the arclength constraint

The Euler-Lagrange operator is applied to the integrand

for which the transversality condition

= 0

becomes

= = 0

because and .

The extremal can be obtained from

the general solution of the Euler-Lagrange equation solved in the case of fixed endpoints. However, the endpoint conditions are now

= 0

where for convenience, we use instead of . These conditions give the equations

and

Subtracting the second from the first, we find , or . Clearly, cannot be a solution because that would make the initial and terminal points coincide. Hence, we have the relationship .

The transversality condition requires that . We interpret this to mean that has a vertical tangent at . This can be implemented analytically by requiring to vanish at . Thus, by implicit differentiation, we ahve

= 0

which, at the point becomes

because . Hence, , and the extremal is given implicitly by

As soon as the transversality condition determined the extremal had a vertical tangent at , we knew it was a semicircle. This equation tells us the center for the circle is and the radius is . We can find the value of from the arc length condition since we know that or . Hence, , and the extremal is

Let us now implement these calculations in Maple. Thus, establish via

> depends(y,x);

Warning, alias or macro ddy defined in terms of y

Warning, alias or macro dy defined in terms of y

Warning, alias or macro y'' defined in terms of y

Warning, alias or macro y' defined in terms of y

>

and recall that is

> F;

>

and the Euler-Lagrange equation, and its first integral, is

> q := Euler_Lagrange(F,x,y);

$q := {y+m*sqrt(1+`y'`^2)-`y'`^2*m/(sqrt(1+`y'`^2)) ...$

>

Therefore, the first integral can be put into the form

> Q1;

>

which is the same differential equation we had to solve in the previous example. Since we know the solution to that equation, we remove the aliases on by

> free(y,x);

>

and recall that the general solution of the Euler-Lagrange equation is

> q8;

>

However, the endpoint conditions are now and , where for convenience, we use instead of . These conditions give the equations

> eq1 := subs(x=0,y=0,q8);
eq2 := subs(x=b,y=0,q8);

>

Subtracting the second from the first, we find

> eq3 := simplify(eq1-eq2);

>

whose solution is

> solve(eq3,b);

>

Clearly, cannot be a solution because that would make the initial and terminal points coincide. Hence, we have the relationship .

> q24 := implicitdiff(q8,x,y);

>

Evaluating the derivative at the point , and setting it to zero, gives

> q25 := subs(x=b,y=0,q24)=0;

>

Recalling that , we have

> q26 := subs(b=2*d,q25);

>

so we conclude that , and hence, that the extremal is

> subs(eq1,d=b/2,k=0,q8);

>

To determine analytically, just from the data of the problem, requires evaluating the arclength integral. Once again, we need to find by implicit differentiation, which we do via