Example 46.1

Consider the function

implemented in Maple as

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

>

where the symbol  is used to represent the derivative ' upon input, and `y'` is used to display the derivative upon output.

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.

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

>

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 that 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 .  Writing these points as

>	P0 := [0,0]; P1 := [1/2,q1]; P2 := [1,exp(1)];

>

we can obtain a representation of the approximating spline as

>	YY := Spline([P||(0..2)],x,degree=1);

>

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

>	plot([YY,Y],x=0..1,color=[red,black], labels=[x,y], labelfont=[TIMES,ITALIC,12], tickmarks=[2,3], title=`Figure 46.5`);

>

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))); `` = evalf(J1);

>

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 := simplify(int(f(x,YY,diff(YY,x)),x=0..1)); `` = evalf(J2);

>

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); `` = evalf(J3);

>

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 .  Writing these points as

>	P0; P1 := [1/3,eval(y[1],q1)]; P2 := [2/3,eval(y[2],q1)]; P3 := [1,exp(1)];

>

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

>	YY := Spline([P||(0..3)],x,degree=1);

>

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

>	plot([YY,Y],x=0..1,color=[red,black], labels=[x,y], labelfont=[TIMES,ITALIC,12], tickmarks=[2,3], title=`Figure 46.6`);

>

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))); `` = evalf(J4);

>

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

If we evaluate  along the linear spline, we obtain

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

>

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.

>

Example 46.2

To approximate the extremal of Example 46.1, that is, the functional

with endpoint conditions

where , 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 approximation becomes

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

>

In Figure 46.7 a graph of the exact solution

>

Y;

>

along with the approximate solution , shows little difference between the two solutions.  We have used color (red for the approximate solution, and black for the exact) 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=[red,black], thickness=[3,5], labels=[x,y], labelfont=[TIMES,ITALIC,12], tickmarks=[2,3], title=`Figure 46.7`);

>

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.

>

>