![9. The integral I(a) = Int[0..2] (2+sin(10a))x^a sin(a/(2-x)) dx depends on the parameter a. What is the value a in [0,5] for which I(a) achieves its maximum?](siam-9.gif){kind=small}
The following solution was supplied by Brian Medley
Siam Challenge Problem 9
I'm not very proud of this one - essentially I just used the same method as
on problem 1 to "iron out" the sine cycles. Near the maximum, the values of
the integral agree to about 20 sig figs. I dealt with this by throwing a
lot of computer time at the problem, rather than doing anything clever. I
did need to use a change of variable for x near 0, as the gradient
approaches infinity, upsetting Simpson's rule there.
The UBASIC program below estimates I(alpha) for three alphas:
0.78593367434, 0.78593367435 and 0.78593367446 (I will rely on the
comments in the program itself to explain how it works). The results
provide evidence that the middle value yields a greater value for the
integral than the others, so I am submitting 0.78593367435 as my answer to
the problem.
Perhaps I should explain how I decided where to search for the optimum
alpha. The factor of (2 + sin(10*alpha)) suggests that the solution will
be near (4n+1)*Pi/20 for some n. To find the right n, I ran the program for
these values (for alpha greater than one I didn't need the change of
variable near 0). It seemed that the place to look was near Pi/4. To see if
this looked right, I wrote an Excel spreadsheet which shows the graph of
the integrand. The value of alpha is controlled by a call to the timer
function, so by holding down the F9 key, you can see the graph changing
with alpha. Alphas near Pi/4 did, as expected, give graphs that looked like
they might yiled the maximum integral.
To narrow down the search for alpha a bit more, I modified the program to
estimate the derivative of I(alpha) and to use linear interpolation to try
to solve I'(alpha)=0. The numerical differentiation involved subtracting
I-values which were closer together than the possible error in I, so this
work was very dubious, but it did lead to an estimate of about 0.7859336744
for alpha, so that is where I looked for the maximum.
Here is the UBASIC program:
10 point 8 ' sets the precision for UBASIC
20 '
30 ' The next line sets up the variables Err and Converg which determine
40 ' respectively the error used in the Simpon estimation and how close
50 ' successive estimates of the integral have to be before we stop doing
60 ' more cycles of the function
70 '
80 Err=10^(-12):Converg=10^(-12):A$="2"
90 '
100 ' The rest of this first bit just sets up a data file for the results
110 '
120 clr time
130 Pi=pi(1)
140 B$=str(int(log(Err*2)/log(10)))
150 S$=A$+B$+".txt"
160 print S$
170 open S$ for output as #1
180 print #1,Err,Converg
190 '
200 ' the next six lines calculate the integral for three values of
210 ' alpha. This is really all the program does. the rest of the code is
220 ' all about the routine for finding the integral
230 '
240 Alph1=0.78593367434:I1=fnI(Alph1)*(2+sin(10*Alph1))
250 Alph2=0.78593367435:I2=fnI(Alph2)*(2+sin(10*Alph2))
260 Alph3=0.78593367446:I3=fnI(Alph3)*(2+sin(10*Alph3))
270 print using(2,26),Alph1;I1
280 print using(2,26),Alph2;I2
290 print using(2,26),Alph3;I3
300 stop
310 '
320 ' Ok, here comes the routine for calculating fnI(alpha), which is the
330 ' integral from 0 to 2 of x(alpha)*sin(alpha/(2-x)
340 '
350 fnI(A)
360 '
370 ' We first do a special calculation for the integral between 0 and
380 ' the point where A/(2-X) equals 3*Pi/2. This is treated separately
390 ' because x^alpha causes the gradient to approach infinity near 0
400 ' for the alphas we are interested in. Simpson's rule can't cope well
410 ' with this so we use the change of variable u=x^alpha.
420 '
430 N=1
440 S=0
450 S0=-1
460 Sest0=100000
470 X0=0
480 X1=fnRoot(3/2)
490 S=fnSimpsonSub(&fnG(),0,X1^A,Err)/A
500 '
510 'start main loop
520 '
530 repeat
540 X00=X0:X0=X1:N=N+1
550 '
560 'find next max or min of sin(A/(2-x))
570 '
580 X1=fnRoot(N+0.5)
590 '
600 'Integrate f(x) from X0 to X1 i.e. a half cycle from a max to a min
610 'or vice versa. Increase the current S (estimate of integral) by this.
620 'Err/N is used for the error as a (fairly shoddily thought out)
630 'safeguard against lots of small errors adding up to a big one.
640 '
650 S+=fnSimpsonSub(&fnF(),X0,X1,Err/N)
660 'print S
670 if N<3 then goto 770
680 '
690 ' for N>3 we can find the ratio of lengths of the last two half cycles
700 ' and use this to do a weighted average of the last two S's, hoping
710 ' to cancel out errors with opposite signs. This average is "Sest"
720 '
730 L=(X1-X0)/(X0-X00)
740 Sest=(S+L*S0)/(1+L)
750 Sest00=Sest0
760 Sest0=Sest
770 S0=S
780 '
790 ' Finally, keep going until successive Sests are within Converg
800 ' of one another.
810 '
820 until abs(Sest00-Sest0)<Converg
830 print using(1,24),A;(2+sin(10*A))*Sest;:print N,time
840 print #1,A;(2+sin(10*A))*Sest;N,time
850 return(Sest)
860 '
870 '
880 ' The routine below finds f(x) to be integrated in the main prog.
890 '
900 fnF(X)
910 return((X^A)*sin(A/(2-X)))
920 '
930 'The Simpson routine below came in a sample program supplied with
940 'UBASIC so I assume it's OK to use it here.
950 '
960 'Simpson's Rule (Numerical Integration)
970 'from: S.Moriguti Suuti Keisanjutu (Kyouritu Shuppan) p199
980 'take x2-x1 small
990 '
1000 fnSimpsonSub(&fnF(),X1,X2,Err)
1010 local Dy,H,M,N,S0,S1,S2,Y1,Y2
1020 N=2:H=(X2-X1)/N
1030 S0=fnF(X1)+fnF(X2):S1=fnF(X1+H):S2=0
1040 Y1=(S0+4*S1)*H/3
1050 repeat
1060 N*=2:H/=2:S2+=S1:S1=0
1070 for M=1 to N step 2:S1+=fnF(X1+M*H):next
1080 Y2=(S0+4*S1+2*S2)*H/3:Dy=Y2-Y1:Y1=Y2
1090 until absmax(Dy)<Err
1100 return(Y2)
1110 '
1120 ' The routine below for fnRoot(n) returns the solution to A/(2-X)=n*Pi
1130 '
1140 fnRoot(N)
1150 return(2-A/(N*Pi))
1160 print fnSimpsonSub(&fnG(),Aa^(A),B^(A),Err),time1000
1170 '
1180 ' The routine below calculates the g(x) to be integrated
1190 ' for the initial, change-of-variable step
1200 '
1210 fnG(X)
1220 return((X^(1/A))*sin(A/(2-X^(1/A))))
Here is the output for various values, of Err and Converg (in fact I
usually took these to be equal).
The first line of output is Err and Converg, written as UBASIC exact
rationals (i.e. using "//" as the fraction bar).
Then, on each line, we have, alph, estimate of I(alpha), number of
iterations of the main loop before estimates differed by less than Converg,
and the running time of the program at which each estimate was produced
(165 hours in one case!). I used different computers for different runs, so
the times do not depend on Err and Converg as you might expect.
With Err = Converg = 10^(-14):
1//100000000000000 1//100000000000000
0.78593367434 3.0337325864855133417014540127080929094 1720 0:01:05
0.78593367435 3.0337325864855133417071924701037717278 1720 0:02:10
0.78593367446 3.03373258648551334106538499431924605367 1720
0:03:15
With Err = Converg = 10^(-18):
1//1000000000000000000 1//100000000000000000000
0.78593367434 3.03373258648549363299366332191643627388 54374
2:44:09
0.78593367435 3.03373258648549363299940044656104040673 54374
5:28:21
0.78593367446 3.03373258648549363235757831051485450566 54374
8:12:31
With Err = Converg = 10^(-20):
1//100000000000000000000 1//100000000000000000000
0.78593367434 3.03373258648549363255021057600563341448 54374
20:51:36
0.78593367435 3.0337325864854936325559477005368710359 54374
41:41:16
0.78593367436 3.03373258648549363255100405981180261064 54374
62:34:40
With Err = Converg = 10^(-21):
1//1000000000000000000000 1//1000000000000000000000
0.78593367434 3.0337325864854936325344063207475530808646215 96692
55:09:36
0.78593367435 3.033732586485493632540143445279265758424267 96692
110:09:36
0.78593367436 3.0337325864854936325351998045513007722682285 96692
165:13:44