Problem 9:

The integral  depends on the parameter .  What is the value  at which  achieves its maximum?

Solution:

This is rather different from Gaston's solution.  The numerical integration works in Maple 8, but not Maple 7.
We want to be able to evaluate
and  so that we can use Newton's method to find the location of the maximum.  Write  so , , and .

We break the interval up into pieces.

On the interval  we use a series in powers of .  On  we use the transformation  , obtaining an integral on .  On the interval  we use Maple's own numerical integration, and on  we expand the integrand in negative powers of y times sin(y).

 > restart: ti:= time():

First for the interval : J1 is the integral, dJ1 its derivative and d2J1 its second derivative.

 > f:= x^alpha*sin(alpha/(2-x)); S:= convert(series(sin(alpha/(2-x)),x,40),polynom): intp:= proc(t) local c,xp,n;   if has(t,x) then xp,c := selectremove(has,t,x);     if xp = x then n:= 1 else n:= op(2,xp) fi;   else n:= 0; c:= t   fi;   collect(expand(c),[cos,sin])*3^(-n-alpha-1)/(n+alpha+1); end: J1:= map(intp,S):

 > dJ1:= diff(J1,alpha): d2J1:= diff(J1,alpha\$2):

 > length(J1),length(dJ1),length(d2J1);

How big is the last term in dJ1 for a typical alpha?

 > evalf(eval(select(has,dJ1,3^(-40-alpha)),alpha=0.8));

Next the transformation.

 > xy := expand(solve(alpha/(2-x)=y,x)); g:= subs(x=xy,f)*diff(xy,y); g1:= eval(g,y=3*alpha/5); ga:= normal(diff(g,alpha)); ga1:= normal(eval(ga,y=3*alpha/5)); gy:= normal(eval(diff(g,y),y=3*alpha/5)); gaa:= normal(diff(g,alpha\$2));

So we want .    is the integral from  to ,  and  its derivatives.  Maple sometimes seems to have trouble evaluating  to high precision, but we don't need as much accuracy there, so the relative error tolerance  can be increased in this case.

 > dH:=diff(int(h(a,y),y=3*a/5 .. 2*Pi),a); d2H:=diff(int(h(a,y),y=3*a/5 .. 2*Pi),a\$2);

 > J2:= a -> evalf(Int(subs(alpha=a,g),y=3*a/5 .. 2*Pi));

 > dJ2:= a -> evalf(Int(subs(alpha=a,ga),y=3*a/5 .. 2*Pi) - 3/5*subs(alpha=a,g1)); d2J2:= a -> evalf(Int(subs(alpha=a,gaa),y=3*a/5 .. 2*Pi,epsilon=10^(-Digits/2)) - subs(alpha=a,6/5*ga1+9/25*gy));

 > J2(0.8),dJ2(0.8),d2J2(0.8);

For the interval , we will use a series for  in negative powers of , relying on Maple to integrate .  We have  where .

 > J3:=add(2^(2+alpha-n)*(-alpha)^n*pochhammer(alpha-n+3,n-3)/(n-2)!*int(sin(y)/y^n,y=2*Pi..infinity),n = 2 .. 30):

 > dJ3:= diff(J3,alpha): d2J3:= diff(J3,alpha\$2):

 > evalf(eval([J3,dJ3,d2J3],alpha=0.8));

 > t30:= add(2^(2+alpha-n)*(-alpha)^n*pochhammer(alpha-n+3,n-3)/(n-2)!*int(sin(y)/y^n,y=2*Pi..infinity),n = 30..30); evalf(eval([t30,diff(t30,alpha),diff(t30,alpha\$2)],alpha=0.8));

 > Jp:= a -> evalf(eval(J1+J3,alpha=a))+J2(a): Ip:= a -> evalf(2+sin(10*a))*Jp(a):  Jp(0.8);

Here it is plotted (with low precision),  in red and  in green.  The maximum is near .

 > plot([Jp,Ip],0..5);

Here is a Newton iteration procedure.  It takes an initial  value, prints the values of ,  and , and returns the new alpha for the next iteration: .

 > Newt:= proc(a)  local J,dJ,d2J,Ia,dI,d2I;  J:= evalf(eval(J1+J3,alpha=a)+J2(a));  dJ:= evalf(eval(dJ1+dJ3,alpha=a)+dJ2(a));  d2J:= evalf(eval(d2J1+d2J3,alpha=a)+d2J2(a));  Ia:= (2+sin(10*a))*J;  dI:=10*cos(10*a)*J+(2+sin(10*a))*dJ;  d2I:= -100*sin(10*a)*J+20*cos(10*a)*dJ+(2+sin(10*a))*d2J;  print([Ia,dI,d2I]);  a-dI/d2I;  end;

 >

 > a1:=Newt(0.8);

 > Digits:= 36: for i from 2 to 5 do a||i:= Newt(a||(i-1)) od; (time()-ti)*seconds;

This has 33 correct digits.

 >