0
$\begingroup$

I need to program in an algorithm that recursively makes algebraic replacements which leads to an utterly complicated algebraic function of $x$, but whose final result is only needed at fixed order in the Laurent series in $x$. I wish to use this fact to make the algorithm go faster. Here is sample problem that mimics the realistic problem:

\begin{align}f_n(x) &= \sum_{j=1}^4\frac{1}{6+j x}n\left(\frac{1}{x}+j \,f_{n-1}(x)\right); \qquad f_0(x) = {\tt f[0]}\\ F(x) &= \sum_{n=1}^6\frac{1}{2+n x}\left(\frac{n}{x}+f(n)\right) \end{align}

The algorithm is to repeatedly insert the first equation into the second one until all instances of $f_n(x)$ turn into $f_0={\tt f[0]}$. All that is needed in the end is the Laurent expansion out to $\mathcal{O}(x^0)$ (just need the $1/x$ and constant terms). nota bene: The sums are there only to lengthen the expression to make Mathematica work harder.

Here is a naive implementation (with the result on my machine):

(*Implementation 1*)
f[n_] /; n > 0 -> Sum[1/(6 + j x) n (1/x + j f[n - 1]), {j, 1, 4}];
Sum[1/(2 + n x) (n/x + f[n]), {n, 1, 6}] //. %;
Timing[Series[%, {x, 0, 0}]]

Result of implementation 1

Performing the actual recursion in the second line is lightening fast. But the final step Series function is painfully slow (4.109 seconds), since the $x$ in the denominators forces Mathematica to appeal to its analytical routines to expand things like $\frac{1}{1+x}\approx 1 -x$.

Since all that I need is the first two terms of the Laurent series, I went ahead and, by hand, expanded the fractions $\frac{1}{6+jx}\approx(\frac{1}{6}-\frac{j x}{36})$ and $\frac{1}{2+nx}\approx(\frac{1}{2}-\frac{n x}{4})$, with the intention of partially freeing Mathematica of painful analytics.

(*Implementation 2*)
f[n_] /; n > 0 -> Sum[(1/6 - (j x)/36) n (1/x + j f[n - 1]), {j, 1, 4}];
Sum[(1/2 - (n x)/4) (n/x + f[n]), {n, 1, 6}] //. %;
Timing[Series[%, {x, 0, 0}]]

enter image description here

This is 8 times faster; the more I save Mathematica from doing analytic manipulations, the faster it goes. But I need it to go faster, yet. How can I carry out the alogorithm to minimize the extent Mathematica has to do computationally costly analytic work?

Here's my failed attempt: ask Mathematica drop higher order terms at each stage in the recursion by adding an O[x] to the end of the line:

(*Implementation 3 [FAILURE]*)
f[n_] /; n > 0 -> Sum[(1/6 - (j x)/36) n (1/x + j f[n - 1]), {j, 1, 4}];
Timing[Sum[(1/2 - (n x)/4) (n/x + f[n]), {n, 1, 6}] + O[x] //. %]

enter image description here

Very fast, but the wrong answer. It gets the wrong answer because it drops terms like x*f[n] as it thinks f[n] is order O[x]^0 when in reality there is a 1/x part. What can I do to speed this up?

$\endgroup$
3
$\begingroup$

Truncate early and often.

f[0] = f0;
f[n_] /; n > 0 := 
  f[n] = Sum[
    Series[1/(6 + j x) n (1/x + j f[n - 1]), {x, 0, 0}], {j, 1, 4}];

Timing[
 res = Series[
   Sum[Series[1/(2 + n x) (n/x + f[n]), {x, 0, 0}], {n, 1, 6}], {x, 0,
     0}]]

(* Out[140]= {0.012000, SeriesData[x, 0, {
Rational[1016635, 162], Rational[5, 486] (-3351289 + 835701 f0)}, -1, 
  1, 1]} *)
$\endgroup$
  • 1
    $\begingroup$ Note to self: the key is not to use ReplaceRepeated but to assign DownValues $\endgroup$ – QuantumDot Oct 1 '14 at 9:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.