

Solutions to Trefethen's problems from
"Hundreddollar, Hundreddigit challenge",
SIAM News volume 35, number 1.
Keith Briggs 2002 Mar 25
Notes on highprecision arithmetic:
http://keithbriggs.info/xrc.html
www.nersc.gov/~dhbailey/mpdist/mpdist.html
www.shoup.net/ntl/
Problem 1

integrate_0^1 cos(log(x)/x)/x dx
Solution:
0.3233674316777787...
Method:
integrate_0^1 cos(log(x)/x)/x dx
= integrate_0^infty W'(y) cos(y) dy
= sum_{n=1; n odd}^infty integrate_{pi}^{pi} W'(u+n*pi) cos(u) du ...(1)
where W=Lambert's W function: W(x)exp(W(x))=x; W'(x)=W(x)/(x(1+W(x))).
The advantage of this representation is that each integral is a product
of cos with slowly varying function (W(x)~log(x)). Thus, it is very well
approximated by an expansion
integrate_{pi}^{pi} W'(u+n*pi) cos(u) du ~=
a_2/2! W^{(3)}(n*pi) + a_4/4! W^{(5)}(n*pi) + ... ...(2)
where a_k = integrate_{pi}^{pi} x^k cos(x) dx.
Both the a_k and the derivatives W^{(k+1)} of W can be computed from
simple recurrences. Already at n=13, two terms of this expansion
give sufficient accuracy. Thus, we do the first few terms in the sum (1)
by numerical integration, then switch to (2) for the rest of the sum.
We may estimate the truncation error as follows: W^{(3)}(x) is asymptotically
2/(x*log(x))^3 for large x. Bounding this by 2/x^3, we can sum the tail
as being < 4*pi/x^2. From this we estimate that 10^6 terms
are sufficient for 10 dp in the final result.
Doing a more accurate tail estimation:
sum_{k=n}^infty W^{(3)}(k*pi)
~= 2*pi sum_{k=n}^infty 1/(x*log(x))^3
~ 2*pi integral_n^infty 1/(x*log(x))^3 dx
= 2*pi [ 2Ei(2log(n)) + (12log(n))/(2n^2 log(n)^2) ]
This suggests that the tail is:
<8e9 at n=1000
<4e11 at n=10000
<2e13 at n=100000
<2e15 at n=1000000
Problem 2

Circular mirror of radius 1/3 at each point of Z^2.
Photon starts at (1/2,1/10) with velocity (1,0).
What is its position at time 10?
Method:
1. Run a quickanddirty python program (trefethen2.py) using
IEEE double arithmetic to get a sequence of bounces, which we
will afterwards confirm with precise arithmetic.
This program will also output data for feeding to "graph" to make
a postscript figure (trefethen2fig.ps).
The sequence of the centres of circular mirrors computed is
[[1, 0], [1, 1], [0, 2], [ 0, 1], [1, 2], [ 0, 0], [1, 0], [0, 1],
[0, 0], [1, 0], [0, 1], [1, 1], [1, 0], [1, 1], [ 1, 0]];
2. Run a C++ program (trefethen2xr.cc) which will recompute
everything in such a way that the final time and velocity vector
are accurate to any requested precision (e.g. 10 decimals).
If the sequence of bounces above was wrong, this program will crash
with a "sqrt of negative argument" error.
Solution:
x(10) = 0.736292698609618310775737171866
y(10) = 0.669642696363571375194376070966
distance = 0.99526291944335416089031180942
More precisely:
x(10) = 0.736292698609618310775737171866391720451575491031024118015570
32135829096701522929400945112122000980773908774827834875197238
77531369315078617517331539949084939398476219513951166249120103
0045569289350602e0
y(10) = 0.669642696363571375194376070967098466862227372928775421597411
09381282909191783522403449125396431985636370847144665160568810
35768639394200390550524425747382398387223992432483782229188924
2655423072136579e0
distance = 0.9952629194433541608903118094267216210294669227341543498032088
580729861796228306320991749818976188760306062973261905117220872
365706973813827798622298601793010235848291712247521533893738542
1202854314982e0
Problem 3

A=
1/1 1/2 1/4 1/7 1/11 ...
1/3 1/5 1/8 1/12 1/17 ...
1/6 1/9 1/13 1/18 1/24 ...
1/10 1/14 1/19 1/25 1/32 ...
1/15 1/20 1/26 1/33 1/41 ...
....
What is the 2norm of A considered as a operator on l^2?
Solution:
1.27422415274...
Method:
We have 1/A_{ij}=(i+1)(i+2(j+1))/2+j(j1)/2 (indexing from 0)
and therefore must maximize over x_j: sqrt(sum_i (sum_j A_{ij}*x_j)^2)
subject to sum_i x_i^2=1
I truncated these sums at n=4000 and performed the constrained
optimization.
Problem 4

minimize exp(sin(50x))
+sin(60exp(y))
+sin(70sin(x))
+sin(sin(80y))
sin(10(x+y))
+(x^2+y^2)/4
Solution:
3.30686864747523728...
at
x=0.0244030796943785 y=2.10612427155358
Method:
Random search followed by gradient optimization.
Check with gp:
41 100
x=0.0244030796962585; y=0.210612427157718
objective:
exp(sin(50*x))+sin(60*exp(y))+sin(70*sin(x))+sin(sin(80*y))sin(10*(x+y))+(x^2+y^2)/4
gradient:
x/2+50*exp(sin(50*x))*cos(50*x)10*cos(10*(x+y))+70*cos(x)*cos(70*sin(x)
y/2+60*exp(y)*cos(60*exp(y))10*cos(10*(x+y))+80*cos(80*y)*cos(sin(80*y))
Problem 5

Find the cubic minimizing supnorm error for 1/gamma(z) on unit disk.
Solution:
min supnorm error = 0.2143352346...
Minimizing cubic:
0.0055418508 + 1.0197619z + 0.62521197z^2  0.60334332z^3
Method:
1. Since 1/gamma(z) has a Taylor series around 0 with real coefficients,
so will the best approximating polynomial.
2. By the maximum modulus theorem, the maximum approximation error
in the unit disk will occur on the boundary.
3. A rough check shows the max error occurs near t=0.2233, 0.36 and at 1/2,
where z=exp(2*pi*i*t).
4. Expand 1/gamma(z) around the first two of these points to get
accurate values to fit to. (1/gamma(1)=0 for t=1/2)
5. Use a gradient optimizer to get the 4 coefficients of the cubic;
to find the max error near the above two approximate t values use a
Brent onedimensional optimizer in the argument t for each gradient step.
Problem 6 NB: this solution is incorrect!

Particle walks on Z^2.
P(up)=1/4, P(dn)=1/4, P(rt)=1/4+eps, P(lt)=1/4eps.
P(eventual return to origin)=1/2.
What is eps?
Solution:
0.108996971771294...
Method:
Compute probabilities on 4000*4000 lattice, use bisection to find eps.
Problem 7

A is a 20000*20000 integer matrix
A(i,i) = prime(i)
A(i,j) = 1 if ij=1,2,4,8,...,16384
= 0 else
What is (A^1)_{11}?
Solution:
0.72507834626840116746868771925116096886918059447950895787816...
Method:
Solve A*x = [1,0,0,0,...,0]^T; solution is x_1.
Conjugate gradient works well since A is sparse symmetric
positive definite. Neither A nor A^1 need be stored.
Problem 8

Square domain [0,2]*[0,2], 3 edges held at 0, one edge (x=2) at 5.
u_t=u_{xx}+u_{yy}
At what time does u at (1,1) reach 1?
Solution:
0.42401138703368836379743366859326...
Method:
Solving the heat equation with Fourier series, we find that the
temperature u at (1,1) at time t is
(5/4) + 40/pi^2 sum_{m=1; m odd}^{infty} sum_{n=1; n odd}^{infty}
(1)^{(m+n)/2} n/m 1/(m^2+n^2) exp[(m^2+n^2)*pi^2*t/4]
where the 5/4 term is the equilibrium temperature.
Bisection easily finds the required value of t.
Problem 9

argmax_{a in [0,5]} (2+sin(10a)) * integrate_0^2 x^a * sin(a/(2x)) dx
Solution:
a = 0.78593367197930...
Method:
The integral is easier to handle if we rewrite it as
integrate_{1/2}^infinity (21/u)^a/u^2 * sin(a*u) du
I used an adaptive integration algorithm optimized for
oscillatory integrands, together with a goldensearch search.
Problem 10

Rectangle R=[5,5]*[1/2,1/2]
Particle starts Brownian motion at (0,0)
What is probability of first exit through a short side of R?
Solution:
2*1.918793989625613051703566593102...e7
= 0.00000038375879792512261034071331862048...
Method:
Unless I have misunderstood something, this problem is easy!
Solve FokkerPlanck equation on R with appropriate boundary conditions,
p(x,y)=prob(first exit through right edge of R starting at (x,y))
which just gives Laplace's equation p_{xx} + p_{yy} = 0
BC: 0 on left, bottom, top edges; 1 on right edge
Solution by separation of variables is
p(0,0) = sum_{m=1, m odd}^{5nfty} 4/(m*pi)*sinh(5*m*pi)/sinh(10*m*pi)
Already 2 terms give almost 30 d.p.
This is the probability of first exit through the righthand edge;
we double this for the probability of first exit through either short edge.

This website uses no cookies. This page was last modified 20130512 10:17
by .

