![1. What is lim[eps->0] Int[eps...1] cos(log(x)/x)/x dx?](siam-1.gif){kind=small}
The following solution was supplied by Brian Medley
Siam Challenge Problem 1
(1)
I firstly used the change of variable x=e^(-u) to transform the integral to:
integral (0 to infinity) (cos(u.e^u)du)
(Actually, I did integral (0 to infinity) (cos(x.e^x)dx); all x's in the
material below are references to the transformed variable, not the original
one)
(2)
Since the cycles of the graphs get narrower and narrower, I split it into
half-cycles (i.e. sections from a maximum to a minimum or vice versa),
integrating each one separately with Simpson's rule. I used UBASIC to allow
me to use fixed-point arithmetic to plenty of decimal places. One of the
programs that came with UBASIC included a Simpson routine so I lazily used
that instead of writing my own.
To find the ends of each half-cycle, I used Newton-Raphson to solve x.e^x=n.pi
As x increases, the cycles of cos(x.e^x) get more equal-sized as x
increases, so in integrating half-cycles from a maximum to a minimum or
vice versa, the positive and negative bits cancel out more and more, so the
current estimate of the integral changes by relatively small amounts.
(3)
There is still quite a bit of wobble between succesive estimates of the
integral, so I used a weighted average of the last two estimates.
I reasoned that the infinite, omitted, part of the graph (which is the main
part of the error) is effectively turned upside down and distorted slightly
at each stage. I approximated the distortion as a horizontal stretch using
the ratio of the lengths of the last two intervals used in the
integration. That is, I assumed that if S_(n-1) and S_n are the last two
estimates of the integral I, then S_(n-1)=I+e and S_n=I-L*e, where L is the
ratio of intervals (nth interval)/((n-1)th interval)
This leads to the weighted average (S_n + L*S_(n-1))/(1+L)
In fact I don't see why I shouldn't average the last two (or more)
estimates in any way I like: the proof of the pudding is in how quickly the
thing converges. I did some trials with averaging the averages and these
did speed up the convergence some more, but they are not included in the
program below.
(4) The program:
10 point 10
20 input Err ' to be used in NewtonRaphson and Simpson's rule
21 S$="si1"+str(int(log(Err)/log(10)))+".txt"
25 clr time
30 '
40 'initialise stuff
50 '
60 Pi=pi(1)
70 X0=0
80 N=1
90 S=0
100 S0=-1
110 X0=0
120 '
130 'start main loop
140 '
150 repeat
160 X00=X0:X0=X1
170 '
180 'find X1 such that f(X1)=n*pi
190 '
200 X1=fnSOLVE(N)
210 '
220 'Integrate f(x) from X0 to X1 i.e. a half cycle from a max to a min
230 'or vice versa. Increase the current S (estimate of integral) by this.
240 'Err/N is used for the error as a (fairly shoddily thought out)
250 'safeguard against lots of small errors adding up to a big one.
260 '
270 S+=fnSimpsonSub(&fnC(),X0,X1,Err/N)
280 if N<3 then goto 380
290 '
300 ' for N>3 we can find the ratio of lengths of the last two half cycles
310 ' and use this to do a weighted average of the last two S's, hoping
320 ' to cancel out errors with opposite signs. This average is "Sest"
330 '
340 L=(X1-X0)/(X0-X00)
350 Sest=(S+L*S0)/(1+L)
355 if abs(Sest-Sest0)<10^(-12) then open S$ for output as #1:print #1,Sest00,Sest0,Sest,time:close:stop
360 print N,Sest
365 Sest00=Sest0
370 Sest0=Sest
380 S0=S
390 N=N+1
400 until 0=1
410 '
420 'routine to calculate f(x)
430 '
440 fnE(X)
450 return(X*exp(X))
460 fnC(X)
470 return(cos(fnE(X)))
480 '
490 'routine to solve xexpx=n*pi
500 'by less than Err, so returned value will be v. close to root.
510 '
520 fnSOLVE(N)
530 local A,C,E
540 'solve xexpx=npi
550 C=N*Pi
560 A=log(C)
570 repeat
580 E=fnE(A):A=A-(E-C)/(E+exp(A))
590 until abs(E-C)<Err
600 return(A)
610 '
620 'The Simpson routine below came in a sample program supplied with
630 'UBASIC so I assume it's OK to use it here.
640 '
650 'Simpson's Rule (Numerical Integration)
660 'from: S.Moriguti Suuti Keisanjutu (Kyouritu Shuppan) p199
670 'take x2-x1 small
680 '
690 fnSimpsonSub(&fnF(),X1,X2,Err)
700 local Dy,H,M,N,S0,S1,S2,Y1,Y2
710 N=2:H=(X2-X1)/N
720 S0=fnF(X1)+fnF(X2):S1=fnF(X1+H):S2=0
730 Y1=(S0+4*S1)*H/3
740 repeat
750 N*=2:H/=2:S2+=S1:S1=0
760 for M=1 to N step 2:S1+=fnF(X1+M*H):next
770 Y2=(S0+4*S1+2*S2)*H/3:Dy=Y2-Y1:Y1=Y2
780 until absmax(Dy)<Err
790 return(Y2)
(5) Results: The last three estimates are given and the time taken for the
program to reach consecutive estimates differing by less than 10^-12, for
various values of Err. In each of these cases, the program terminated at
N=4448
Err=10^-10
0.323367431677898988394759195828382377394353978578
0.323367431678898989038286816431033685175460732908
0.32336743167789965874451240775196565083512804012 0:00:57
Err=10^-14
0.323367431677278678606420806987337147771805233512
0.323367431678278679153090031163620197210020231264
0.323367431677279348956109089947779145402335922849 0:08:10
Err=10^-18
0.323367431677278593457282291807802601493038273706
0.323367431678278594003950039787285108249009474688
0.323367431677279263806970573778677032689285130044 1:13:31