Unit 9: Calculus of Variations

Chapter 46: Basic Formalisms

Section 46.2: direct methods

Copyright

Copyright * 2001 by Addison Wesley Longman, Inc.

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;

> with(plots):

Warning, the name changecoords has been redefined

>

Introduction

In this section we present two methods for approximating , the extremal which satisfies the endpoint conditions

and which stationarizes the functional

First, we examine a finite difference technique of Euler, and then, we consider the Rayleigh-Ritz procedure. Such methods are called direct methods , as opposed to the indirect method of forming and solving the Euler-Lagrange differential equation, which will be seen in the next section.

>

Euler's Method of Finite Differences

Euler's method of finite differences begins by discretizing the interval with equispaced nodes

where . The function-values are designated as the unknowns . The derivatives y' are replaced with the quotients

The integral is then approximated by the finite sum

which becomes a function of the unknowns . Of course, = , and = .

The variables are then determined by the techniques of calculus.

>

Example 46.1

Consider the function

implemented in Maple as

> f := (x,y,yp) -> y^2+yp^2+2*y*exp(x);

>

where the symbol is used to represent the derivative '.

Let the functional be given by

with the endpoint conditions being

The exact solution to this calculus of variations problem is the function

> Y := 1/2*exp(1)*csch(1)*sinh(x)+x/2*exp(x);

>

a solution we will learn to calculate in Sections 46.3 and 46.4. (See Example 46.7.)

The graph of the exact solution is given by

> plot(Y,x=0..1, color=black);

>

That this solution satisfies the endpoint conditions is seen from the following calculations.

> simplify(subs(x=0,Y));
simplify(subs(x=1,Y));

>

Since = 1, we have , and . For the moment, is left unspecified to allow for increasingly more accurate solutions as is increased.

The function is given by

> phi := n -> sum(f(k/n,y[k],(y[k+1]-y[k])/(1/n))*(1/n),k=0..n-1);

>

where we have made it a function of the index . Thus, for example, if , the sum which approximates the functional is

> phi(2);

>

Notice that the sum includes the endpoint values and = . These endpoint values are know, so the approximating sum becomes

> q := subs(y[0]=0,y[2]=exp(1),phi(2));

>

There is but the one unknown, namely, , corresponding to . It is determined by setting to zero the derivative with respect to . This yields

> q1 := solve(diff(q,y[1]),y[1]);

>

The extremal is being approximated by a linear spline connecting the three points . If the and coordinates of these points are entered into the lists

> Lx := [0,1/2,1];
Ly := [0,q1,exp(1)];

>

we can obtain a representation of the approximating spline as

> YY := spline(Lx,Ly,x,linear);

>

A graph of the approximate extremal, in cyan, is compared to the exact extremal, in black, in Figure 46.5.

> plot([YY,Y],x=0..1,color=[cyan,black], labels=[x,y], labelfont=[TIMES,ITALIC,12], xtickmarks=2, ytickmarks=3);

>

There are two ways the value of the functional might be approximated. First, its value can be taken as the value of the discrete sum , obtained in Maple with

> J1 := simplify(subs(y[0]=0,y[1]=q1,y[2]=exp(1),phi(2)));

>

Since the known values of and appear in the sum , it is necessary to supply those values when evaluating the sum.

An alternate approximation of can be computed by evaluating it along the piecewise linear path just computed as the approximating linear spline. In this event, we would compute

> J2 := int(f(x,YY,diff(YY,x)),x=0..1);

>

We can also obtain the value of along the exact extremal. This turns out to be

> J3 := collect(int(f(x,Y,diff(Y,x)),x=0..1),csch);

>

The floating-point equivalents for all three values of are

> evalf(J1);
evalf(J2);
evalf(J3);

>

The first is least accurate. That is because in the passage from to , the extremal is approximated by a degree zero spline. The value reflects this use of a lower-degree approximation of the extremal. After the values of the approximating extremal are obtained at the nodes, a linear spline is formed. Hence, will clearly be a better approximation to the exact stationary value given by because is computed by integrating along the more accurate linear spline, and not the degree zero spline used to obtain .

If we next take , we will have two interior nodes at which to compute the solution. The approximating sum for the case is

> q := subs(y[0]=0,y[3]=exp(1),phi(3));

>

where we have supplied the known values and . The two unknowns and are determined by the ordinary techniques of multivariable calculus. Thus, set to zero the derivatives with respect to each of the unknowns, and solve, obtaining'

> q1 := solve({diff(q,y[1]),diff(q,y[2])},{y[1],y[2]});

>

The extremal is being approximated by a linear spline connecting the four points . If the and coordinates of these points are entered into the lists

> Lx := [0,1/3,2/3,1];
Ly := [0,op(subs(q1,[y[1],y[2]])),exp(1)];

>

we can obtain the following representation of the approximating linear spline.

> YY := spline(Lx,Ly,x,linear);

>

A comparison between this more accurate spline and the exact extremal is given in Figure 46.6, where the spline is drawn in cyan, and the exact extremal, in black.

> plot([YY,Y],x=0..1,color=[cyan,black], labels=[x,y], labelfont=[TIMES,ITALIC,12], xtickmarks=2, ytickmarks=3);

>

Approximating the value of the functional with the discrete sum gives

> J4 := simplify(subs(y[0]=0,op(q1),y[3]=exp(1),phi(3)));

>

where again, the values and are needed in this sum.

If we evaluate along the linear spline, we obtain

> J5 := int(f(x,YY,diff(YY,x)),x=0..1);

>

Once again, we look at the floating-point equivalents

> evalf(J4);
evalf(J5);
evalf(J3);

>

to see that evaluating the functional along the linear spline gives a more accurate value than , obtained by evaluating the functional along a degree-zero spline.

>

The Rayleigh-Ritz Technique

The Rayleigh-Ritz technique for approximating extremals in the calculus of variations bears some resemblance to the Rayleigh-Ritz approximation scheme (Section 15.3) for solving differential equations. Texts such as [1] point out that the name is sometimes shortened to just the Ritz method since W. Ritz greatly generalized, in 1908 and 1909, the earlier work of Lord Rayleigh.

To approximate the extremal which stationarizes the functional

with respect to curves satisfying the conditions

replace with the sum + , where satisfies the nonhomogeneous boundary conditions

and = 0, . This makes become a function of the parameters , and the ordinary techniques of multivariable calculus are used to determine values which starionarize . (Note that the same was used in Section 15.3.)

>

Example 46.2

To approximate the extremal which stationarizes the functional

with endpoint conditions

where , that is, the functional of the previous section, we take as

> Phi := x*exp(x)+add(c[k]*x*(x^k-1),k=1..2);

>

The function satisfies the nonhomogeneous boundary conditions, and the functions , satisfy the homogeneous conditions = 0, . This converts the functional to

> Jc := Int(f(x,Phi,diff(Phi,x)),x=0..1);

>

with value

> Q := value(Jc);

>

To determine the parameters and , set to zero the derivatives with respect to each parameter, and solve, obtaining

> q := solve({diff(Q,c[1])=0,diff(Q,c[2])=0},{c[1],c[2]});

>

Hence, the Rayleigh-Ritz approximant becomes

> Phi[1] := subs(q,Phi);

>

A graph of the exact solution

> Y;

>

along with the approximate solution , shows little difference between the two solutions. We have used color (cyan for the approximate solution, and black for the exact) and varied the line thickness (thick for ) to distinguish one curve from the other. However, they are so close as to appear as virtually the same curve.

> plot([Phi[1],Y],x=0..1, color=[cyan,black], thickness=[3,1], labels=[x,y], labelfont=[TIMES,ITALIC,12], xtickmarks=2, ytickmarks=3);

>

The values of on each curve are computed by

> Je := evalf(Int(f(x,Y,diff(Y,x)),x=0..1));
Jc := evalf(Int(f(x,Phi[1],diff(Phi[1],x)),x=0..1));

>

The functional , computed on the approximate extremal, is marginally larger than the value on the exact extremal.

>

References

[1] Marvin J. Forray, Variational Calculus in Science and Engineering, McGraw-Hill Book Company, 1968.

>