Home > Uncategorized > SIAM 100-digit challenge no 1 – Matlab and Pari GP

SIAM 100-digit challenge no 1 – Matlab and Pari GP

Continuing the SIAM 100-digit challenge presentation I show below how to solve the first problem of the series. A great bibliography for the series is the book entitled The SIAM 100-digit Challenge. A study in high-accuracy numerical computing by Folkmar Bornemann, Dirk Laurie, Stan Wagon, Jorg Waldvogel. The first problem in the series deals with the numerical integration of a highly oscillating function. As stated, the purpose is to obtain at least ten significant digits in the numerical computation.

Problem 1. Compute the limit

$\displaystyle \lim_{\varepsilon \rightarrow 0} \int_\varepsilon^1 x^{-1}\cos(x^{-1}\log(x))dx.$

The first step is to plot the function which is being integrated in order to have an idea of what we are dealing with.

Seeing these oscillations near ${0}$ is not too reassuring… There are numerical techniques which can deal with these highly oscillatory integrals. One is to split the interval into many smaller intervals accumulating at ${0}$ such that on each interval the function is simple enough such that a simple quadrature rule can be used. I will not be doing this.

Instead, let’s notice that the function does not have any singularities, except in zero, if we allow ${x}$ to take also complex values. Thus, Cauchy’s theorem regarding complex intervals says that the value of the integral does not change when we take different paths from the starting point to the end point. If ${z :[0,1] \rightarrow [0,1]}$ is defined such that ${z(0)=0}$ and ${z(1)=1}$ then the integral is equal to

$\displaystyle \int_0^1 \text{Re}(z(t)^{i/z(t)-1}z'(t))dt.$

In the experiment below I took ${z(t)=t+i\sin(\pi t)}$. First, let’s notice that the graph of this new function is much more nicer and we expect that standard integration techniques would give satisfying results.

Here’s a Matlab code that can find the first ${15}$ digits correctly. On the site of Nick Trefethen you can find the answers up to ${40}$ digits and you can check that this result has the correct first ${15}$ digits.

function problem1siam
x = 0.01:0.0001:1;
figure(1)
g = @(t) cos(log(t)./t)./t;
y = g(x);
plot(x,y)
figure(2)
z = h(x);
plot(x,z);
function res = h(t)
z = t+i*sin(pi*t);
res = real(z.^(i./z-1).*(1+i*cos(pi*t)*pi));


The Matlab result is

 0.323367431677779


The question is now, can we do more? Unfortunately Matlab cannot do multiprecision computations by default. All Matlab result are given to machine precision, that is 16 significant digits. Additional non-free toolboxes need to be installed. A nice free software which can do great multiprecision computations is Pari GP. It turns out that it has a numerical integration routine that we can put to the test below. We set the precision to ${100}$ digits and we find a result. All digits shown are correct, as we can see by comparing with the results in the site linked in the beginning. At least the first ${40}$ digits correspond to the ones on Nick Trefethen’s site.

\p 100;

g(t)    = if(t==0,return(0),\
z = t+I*sin(t*Pi); return(real(z^(I/z-1)*(1+I*cos(t*Pi)*Pi))));

Result = intnum(x=0.0000000000000001,1,g(x));


The 100 digit result is the following:

0.3233674316777787613993700879521704466510466257254696616810364434317903372106728944319303704641024514