Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

I want it to be stable near $f(0) = 1$. Is there a nice function that does this already, like maybe a hyperbolic trig function or something like expm1, or should I just check if $x$ is near zero and then use a polynomial approximation?

share|improve this question
1  
Once you use expm1 to compute $\exp(x)-1$, there's no further loss of significance in dividing $x$ by it, is there? –  Rahul Aug 21 '12 at 22:13
1  
@Rahul: You are suggesting to use x/expm1(x) but with a single extra check for 0/0 when x is exactly 0? –  martin Aug 21 '12 at 22:29
    
Suggesting, yes, (though to be clear, not with a great deal of authority). –  Rahul Aug 21 '12 at 23:06
1  
Why not Bernoulli numbers? –  Pedro Tamaroff Aug 22 '12 at 0:46
    
$$f(x)={e^{-x/2}\over{\rm sinc}(i x/2)}\ .$$ –  Christian Blatter Aug 22 '12 at 10:09

4 Answers 4

If you don't want to use the expm1() function for some reason, one possibility, detailed in Higham's book, is to let $u=\exp\,x$ and then compute $\log\,u/(u-1)$. The trick is attributed to Velvel Kahan.

share|improve this answer
    
I am totally OK with using expm1(). –  martin Aug 21 '12 at 23:44
    
Does this give any better precision? $u-1=\exp(x)-1$ would seem to lose precision near $x=0$. The computation of $\exp(x)$ is the place where precision is lost. Am I missing something? –  robjohn Aug 22 '12 at 0:54
1  
@robjohn: Neither $u-1$ nor $\log u$ is very accurate for $ u \sim 1$. But if you follow the link, it asserts the errors compensate for each other. –  Hurkyl Aug 22 '12 at 1:54
    
@rob, yes, there is an example there in Higham's book that shows how the effect of the errors in computing the numerator and denominator cancel out. I'll try to find the original ref by Kahan, and will link to it in the answer if I manage to do find the ref. –  J. M. Aug 22 '12 at 11:36
    
@J.M.: Yes, I believe this shows it: $$\frac{\log(\exp(x)+\delta)}{\exp(x)+\delta-1} -\frac{\log(\exp(x)}{\exp(x)-1}\sim\frac12\delta$$ for small $x$ and $\delta$. Very nice. –  robjohn Aug 22 '12 at 12:24

Consider the Bernoulli numbers, defined by the recursive formula:

$$B_0=1$$

$$\sum_{k<n} {n\choose k }B_k=0\text{ ; } n\geq 2$$

This gives the sequence:

$$\{B_n\}_{n\in \Bbb N}=\left\{ 1,-\frac 1 2,\frac 1 6 ,0,\frac 1 {30},0,\dots\right\}$$

It's generating function is

$$ \sum_{n=0}^\infty B_n \frac{x^n}{n!}=\frac{x}{e^x-1}$$

It's first few terms are

$$1-\frac x 2 +\frac {x^2}{12}-\frac{x^4}{720}+\cdots$$

The numbers' denominators grow pretty fast, so you should have no problem with convergence: in fact, the function is essentialy $=-x$ for large negative $x$ and $=0$ for large positive $x$, so a few terms should suffice the "near origin" approximation.

share|improve this answer
    
This is pretty cool and in fact it's what I meant by "check if x is near zero and then use a polynomial approximation." –  martin Aug 22 '12 at 1:57
    
@kxx, you might even consider building a Padé approximant from the series expansion Peter has shown you. For instance, the series expansion of $$\frac{(x-4)^2+4}{\frac{(x+3)^2}{3}+17}$$ matches the original function up to the fourth series term. –  J. M. Aug 22 '12 at 11:34

You mention using hyperbolic functions. You might try $$ \frac{x}{\exp(x)-1}=\frac{x/2}{\exp(x/2)\sinh(x/2)} $$ This loses no precision if the $\sinh$ is computed to full precision by the underlying system.

Note that $\mathrm{expm1}(x)=2\exp(x/2)\sinh(x/2)$.

Example: 15 digit calculations

$x=.00001415926535897932$: $$ \begin{align} \frac{x}{\exp(x)-1} &=\frac{.00001415926535897932}{1.000014159365602-1}\\ &=\color{#C00000}{.9999929203}73447 \end{align} $$ $$ \begin{align} \frac{x/2}{\exp(x/2)\sinh(x/2)} &=\frac{.00000707963267948966}{1.000005000012500\cdot.00000707963267954879}\\ &=\color{#C00000}{.999992920384028} \end{align} $$

share|improve this answer

If your system provides expm1(x) they should have worried about errors and stability at least as well as you can. It is a better question if that is not available. Wolfram alpha gives $1-\frac {x}2+\frac {x^2}{12}$ for the second order series around $0$, so you could check if $x$ is close to zero and use that.

share|improve this answer
    
My system provides expm1. So you are seconding Rahul's suggestion to use x/expm1(x) with the check for x exactly equal to zero? –  martin Aug 21 '12 at 23:46
    
@kxx: that's what I would do. You could try some cases, comparing against Alpha (which will give arbitrary precision), for confirmation. –  Ross Millikan Aug 21 '12 at 23:57

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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