Keith Briggs   ```Solutions to Trefethen's problems from "Hundred-dollar, Hundred-digit challenge", SIAM News volume 35, number 1. Keith Briggs 2002 Mar 25 Notes on high-precision 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)) + (1-2log(n))/(2n^2 log(n)^2) ] This suggests that the tail is: <8e-9 at n=1000 <4e-11 at n=10000 <2e-13 at n=100000 <2e-15 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 quick-and-dirty 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 (trefethen2-fig.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 (trefethen2-xr.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 2-norm 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(j-1)/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 sup-norm error for 1/gamma(z) on unit disk. Solution: min sup-norm 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 one-dimensional 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/4-eps. 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 |i-j|=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/(2-x)) dx Solution: a = 0.78593367197930... Method: The integral is easier to handle if we rewrite it as integrate_{1/2}^infinity (2-1/u)^a/u^2 * sin(a*u) du I used an adaptive integration algorithm optimized for oscillatory integrands, together with a golden-search 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...e-7 = 0.00000038375879792512261034071331862048... Method: Unless I have misunderstood something, this problem is easy! Solve Fokker-Planck 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 right-hand edge; we double this for the probability of first exit through either short edge. ``` This website uses no cookies. This page was last modified 2013-05-12 10:17 by .