In Kahan's account of the Muller's recurrence (p.16): $x_{n+1} = E(x_n, x_{n-1})$ for the function $$ E(y, z) = 108 - \frac{815-\frac{1500}{z}}{y}, $$ he uses the characteristic polynomial for the recurrence to deduce the closed form: $$ x_n = \frac{\alpha 3^{n+1} + \beta 5^{n+1} + \gamma 100^{n+1}}{\alpha 3^{n} + \beta 5^{n} + \gamma 100^{n}}. $$ If $x_0=4, x_1=4.25$ exactly, then we can choose $\alpha=\beta=1, \gamma=0$. But, if there is some numerical round-off in a numerical calculation of the $x_n$, then the closed form becomes ("at least initially"), $$ x_n = \frac{3^{n+1} + 5^{n+1} + \gamma_n 100^{n+1}}{3^{n} + 5^{n} + \gamma_n 100^{n}}\\ = 100 - \frac{95+97(3/5)^n}{20^n\gamma_n+1+(3/5)^n} $$ where $\gamma_n$ is small.

I can't derive the last equality, and can't see how this changes the limit (as $n\rightarrow \infty$) to 100 instead of 5. Furthermore, what justifies the choice of $\alpha$, $\beta$, $\gamma$ above amongst the presumably infinite options?


$\alpha = 1, \beta = 1, \gamma = 0$ gives the exact solution from the starting values. For the round-off error setting they are slightly disturbed; we still have $\alpha \approx 1, \beta \approx 1,$ but $\gamma$ has the order of the machine epsilon $\gamma \approx \epsilon.$ To compute Kahan's form from

$$x_n = \frac{3^{n+1} + 5^{n+1} + \gamma_n 100^{n+1}}{3^{n} + 5^{n} + \gamma_n 100^{n}}$$ first do some algebraic manipulations with the numerator $$3^{n+1} + 5^{n+1} + \gamma_n 100^{n+1}= 3\cdot 3^{n} + 5\cdot 5^{n} + 100\gamma_n 100^{n}$$ $$=(100-97)\cdot 3^{n} + (100-95)\cdot 5^{n} + 100\gamma_n 100^{n}$$ $$=100(3^{n} + 5^{n} + \gamma_n 100^{n})-97\cdot 3^{n} -95\cdot 5^{n} $$ hence the intermediate value for $x_n$ is $$x_n = 100 - \frac{97\cdot 3^n + 95\cdot 5^n}{3^n + 5^n + \gamma_n 100^n}$$ now divide by $5^n$ and get $$x_n = 100 - \frac{ 95+ 97\cdot (3/5)^n }{ 20^n \gamma_n+ 1 +(3/5)^n }$$ Now it is clear that the limit is $\lim_\limits{n\rightarrow \infty} x_n = 100$ for $\gamma_n \ne 0$ and $100-95 = 5$ for $\gamma_n = 0$

The main difference between exact (error=0) and round-off is the growing error term $20^n \gamma_n$. For double precision calculations with a constant $\gamma_n = \epsilon = 2.22\times 10^{-16}$ you have $20^{15}\epsilon \approx 7275.9576$ and $20^{20}\epsilon \approx 23283064365.4$

  • $\begingroup$ Thank you @gammatester! I understand the reasoning much more clearly now. It may be a separate question, but how can we be sure that $\alpha=\beta=1, \gamma=0$ is the correct and only choice? How can this argument be made if a different choice is made? $\endgroup$
    – lauren96
    Nov 11 '17 at 8:42
  • 1
    $\begingroup$ $\alpha$ and $\beta$ cannot both be zero, because this would contradict $x_0=4$. So WLOG you can assume $\alpha=1$ and get two linear equations for $\beta, \gamma$ with the solution $\beta=1, \gamma=0.$ $\endgroup$ Nov 11 '17 at 10:28

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.