## New Algorithm to Make Short Work of Challenging Problems

October 1, 2010 — Jon McLoone, International Business & Strategic Development

Buried deep in the list of new technology in the *Mathematica* development pipeline was the item “integration of oscillatory functions (univariate, multivariate)—new algorithm”. I expect most people will overlook it, as I did, in favor of the new functions, new directions, big infrastructure, and the eye candy. Even worse, most people who will use it won’t even know—it will be selected automatically when needed, like many of *Mathematica*‘s algorithms. So I think it’s my duty to share my discovery that this algorithm is actually really cool.

Why is it so cool?

The first clue I had was when I read in the notes that this was the first time anyone had fully automated the algorithm into a very wide class of problems. Second, that it was a hybrid numeric-symbolic method (putting it beyond the reach of most numerical systems). And finally, that it was developed by the talented Wolfram Research developer Andrew Moylan.

I decided to look for an example to test it on and remembered the SIAM hundred-dollar, hundred-digit challenge. This was a set of “tricky” numerical problems posed by Nick Trefethen of Oxford University, and prize money was awarded for correct solutions. Ninety-four teams from around the world entered, but only 20 of them managed to solve all the tasks.

The first task could be easily stated in *Mathematica* notation:

But *Mathematica* 7, like other systems that do numerical integration, fails on this tricky problem due to the strong oscillations near the origin.

People who solved this challenge had to be much more ingenious and transform the problem through a change of variables into a numerically more stable form, or into a contour integral, and then work on computing with those. It was, after all, meant to be a challenge.

But with this new algorithm, it just works:

A world-class problem that stumped many top minds that you can solve just by typing the question—now that’s cool.

## 9 Comments

Dear Dr.McLoone,

Can you please let me know your email address?

Regards

Ghorbani

I agree this is cool. What would interest me is how it is done but I assume this is meant to be a trade secret. Still nice preview snippet. Can’t wait to get 8.0 into my hands.

One thing we did for these problems was use interval methods to PROVE that the digits were correct. We were able to do that for 7 of the 10 problems…. but not for this one — Problem #1. Basically, we didn’t consider a problem really solved until we could A. get 10000 digits (which we could do for 9 of them) and B. prove that at least the 10 digits were correct (which we could do for, if I recall, 7 of them).

Good progress – although, one must admit that the problems that become “difficult” are always a bit contrived.

I tried to calculate it on Mathematica 7.0.1. It’s kind of hard. My strategy was to define a substitution, x=x(y), such that for small x or y, the argument of the sine differs approximately by pi.

Then the integral of f(x) from 0 to 1 is the same as the integral of f(x(y)) dx/dy over y from 1 to maxy, which is the inverse function of x(y) at argument x=1. However, the two integrands, when identifying the names of the variables, also tend to cancel each other for very small x or y.

However, Mathematica 7.0.1 seems to have a problem to integrate the integrand even though the leading cos/x sub-hyperbolic sine is cancelled. The fast oscillation itself seems to be a problem for it.

I would probably have to go further, to second order, and cancel the function so that only x*cos is left… And I won’t give it that much time. ;-)

Andrew Moylan is at it again. I remember his thought-provoking post on simulating the world cup knock out stage. That was a work of genius.

So Mathematica 7 is able to solve Nick’s integral, you simply need the right approach! Regards Robert Reynolds

I have purchased the Radeon HD 5970 GPU Graphics card and I am anxiously awaiting Mathematica 8′s OpenGL and CUDA capabilities to use with my card.

Again using Mathematica 7, Laurie’s representation and my method we get

laurie[z_] := z^(I/z – 1)

-Re[With[{r = 19.36399},

NIntegrate[

laurie[r Exp[I \[Theta]]] I r Exp[I \[Theta]], {\[Theta], 0,

2 \[Pi]}]]]

0.323367

Wow that really is cool