Computing functions with discontinuity

Question

what is the best way to compute the value of the following function in a computer program, e.g., C/Fortran, in double precision?

$$ g(a,b)={{e^a/a - e^b/b}\over{a-b}} + {1\over{ab}}. $$

The above function has removable discontinuities at a=0, b=0, and a=b. I also need to code a similar three-parameter function that has removable discontinuities at a=0, b=0, c=0, a=b, b=c, c=a, and a=b=c. I am trying to avoid too many if-else in the program, and wondering if there are any mathematical tricks that can help.

As for the one-parameter function $g(x)=(e^x-1)/x$, it is currently coded as follows:

If |x| > $\epsilon$, return g(x) directly; else return $e^{x/2}{{sinh(x/2)}\over{(x/2)}}$, where $sinh(t)\over{t}$ is computed from a truncated Tayler series.

This method can easily be applied to function $f(a,b)=(e^a-e^b)/(a-b)$, as it is equal to $e^{(a+b)/2}{{sinh((a-b)/2)}\over{(a-b)/2}}$, but the best solution for $g(a,b)$ is still not obvious to me.

Thanks.

Answer 1

The answer is the uppermost rightmost element of the matrix exponential $$ \exp\begin{pmatrix} 0 & 1 & 0 \\ 0 & a & 1 \\ 0 & 0 & b \end{pmatrix} $$ So you can use the Taylor series of $e^x$ and apply it to $$ A= \begin{pmatrix} 0 & 1 & 0 \\ 0 & a & 1 \\ 0 & 0 & b \end{pmatrix} $$ Proof: We have $$ A= SDS^{-1} = \\ \phantom{x} \\ \frac{1}{ab(a-b)} \begin{pmatrix} 1 & 1 & 1 \\ 0 & a & b \\ 0 & 0 & -b(a-b) \end{pmatrix} \begin{pmatrix} 0 & 0 & 0 \\ 0 & a & 0 \\ 0 & 0 & b \end{pmatrix} \begin{pmatrix} ab(a-b) & -b(a-b) & a-b \\ 0 & b(a-b) & b \\ 0 & 0 & -a \end{pmatrix} $$ Therefore $$ \exp(A) = S \begin{pmatrix} e^0 & 0 & 0 \\ 0 & e^a & 0 \\ 0 & 0 & e^b \end{pmatrix} S^{-1} $$ and the uppermost rightmost element of this is $$ \frac{1}{ab(a-b)} \begin{pmatrix} 1 & 1 & 1 \end{pmatrix} \begin{pmatrix} e^0 & 0 & 0 \\ 0 & e^a & 0 \\ 0 & 0 & e^b \end{pmatrix} \begin{pmatrix} a-b \\ b \\ -a \end{pmatrix} $$ which is $$ \frac{(a-b)e^0 +be^a-ae^b}{ab(a-b)} = \frac{e^a/a-e^b/b}{a-b} +\frac{1}{ab} $$ Note that the evaluation of the Taylor series can be simplified due to the zeroes in $A$.

We can show that the same Taylor series can also be used in the special cases $a=0$, $b=0$ and $a=b$. So we do not have any distinction of cases at all.

Edit

While the approach described above avoids all distinctions of cases, it might not be the best one for all cases. If your main goal is numeric stability, you should at least distinguish $5$ cases, as described below.

In the following, I use $\mathrm{sinhc}(x)=\frac{\sinh x}{x}$

1. $|a|>\epsilon$, $|b|>\epsilon$ and $|a-b|>\epsilon$. In this case, simply use the original definition of $g(a,b).$
2. $|a-b|\leq \epsilon$, $|a|>\epsilon$ and $|b|>\epsilon$. Rewrite $g$ as follows \begin{align} g(a,b) & = \frac{be^a-ae^b}{ab(a-b)}+\frac{1}{ab} \\ & = \frac{be^a-ae^b + be^b - be^b}{ab(a-b)}+\frac{1}{ab} \\ & = \frac{b(e^a-e^b) + e^b(b-a)}{ab(a-b)}+\frac{1}{ab} \\ & = \frac{1}{a} \cdot \left(\frac{e^a-e^b}{a-b} - \frac{e^b-1}{b}\right) \\ & = \frac{1}{a} \cdot \left(\exp\left(\frac{a+b}{2}\right)\mathrm{sinhc}\left(\frac{a-b}{2}\right) - \frac{e^b-1}{b}\right) \\ & = \frac{1}{b} \cdot \left(\exp\left(\frac{a+b}{2}\right)\mathrm{sinhc}\left(\frac{a-b}{2}\right) - \frac{e^a-1}{a}\right) \end{align}
3. $|a-b|> \epsilon$, $|a|\leq\epsilon$ and $|b|>\epsilon$. Rewrite $g$ as follows: \begin{align} g(a,b) & = \frac{be^a-ae^b}{ab(a-b)}+\frac{1}{ab} \\ & = \frac{be^a-ae^b+a-b}{ab(a-b)} \\ & = \frac{b(e^a-1) - a(e^b-1)}{ab(a-b)} \\ & = \frac{1}{a-b}\cdot\left(\frac{e^a-1}{a} - \frac{e^b-1}{b} \right) \\ & = \frac{1}{a-b}\cdot\left(\exp\left(\frac{a}{2}\right)\mathrm{sinhc}\left(\frac{a}{2}\right) - \frac{e^b-1}{b} \right) \end{align}
4. $|a-b|> \epsilon$, $|a| > \epsilon$ and $|b| \leq \epsilon$. Rewrite $g$ as follows: \begin{align} g(a,b) & = \frac{1}{a-b}\cdot\left(\frac{e^a-1}{a} - \frac{e^b-1}{b} \right) \\ & = \frac{1}{b-a}\cdot\left(\exp\left(\frac{b}{2}\right)\mathrm{sinhc}\left(\frac{b}{2}\right) - \frac{e^a-1}{a} \right) \end{align}
5. None of the above, i.e. at least two of the values $|a|$, $|b|$ and $|a-b|$ are smaller than or equal to $\epsilon$ (which means that the third one is smaller than or equal to $2\epsilon$). In this case, you can use the matrix method proposed above. It will converge quickly due to the small values of $a$ and $b$. However, I do not know how many terms in the truncated Taylor series should be used.

Answer 2

Starting from @Rene Girard's answer, if you want to use Taylor series $$g(b,b)=\frac{e^b (b-1)+1}{b^2}=\sum_{n=0}^\infty \frac{n+1}{(n+2)!} b^n=\sum_{n=0}^p \frac{n+1}{(n+2)!} b^n+\sum_{n=p+1}^\infty \frac{n+1}{(n+2)!} b^n$$

$$S_2=\sum_{n=p+1}^\infty \frac{n+1}{(n+2)!} b^n=\frac{(b+p+2) b^{p+3}+(b-1) e^b (\Gamma (p+4)-\Gamma (p+4,b))}{b^2 \Gamma (p+4)}$$

If $b$ is small we have an underestimate $$S_2 < \frac{(p+2) }{ (p+3)!}b^{p+1}$$ so, at least, we can evaluate $p$ such that $$ \frac{(p+2) }{ (p+3)!}b^{p+1} \leq \epsilon \implies (p+3)! \geq \frac{(p+2)\,b^{p+1}} \epsilon\implies (p+2)! \geq \frac{b^{p+1}} \epsilon$$ This would give $$p \sim b \,e^{1+W(t)}-\frac 52 \qquad \text{where} \qquad t=\frac{1}{2 e b}\log \left(\frac{2 \pi }{b \epsilon ^2}\right)$$ Using $b=2$ and $\epsilon=10^{-10}$, this simple formula gives $p=16.4216$, while $S_2=10^{-10}$ gives exactly $p=15.2549$.

Answer 3

Replacing $a$ by variable $x$ and considering $b$ as a constant, the function g(a,b) can written as follows:

$$ g(x,b) = \frac{1}{xb} \left(\frac{be^x -xe^b}{x-b} + 1 \right) $$

Now consider the limit $x \rightarrow b$ of $g(x,b)$:

$$\lim_{x \to b} \frac{1}{xb} \left(\frac{be^x -xe^b}{x-b} + 1 \right) = \frac{1}{b^2} \left(\frac{be^b -be^b}{b-b} + 1 \right) = \frac{1}{b^2} \left(\frac{0}{0} + 1 \right) $$

Using l'Hospital rule for

$$\lim_{x \to b} \frac{be^x -xe^b}{x-b} = \lim_{x \to b} \frac{be^x-e^b}{1} = e^b(b-1)$$

So

$$g(b,b) = \frac{1}{b^2}(e^b(b-1) +1) $$

Caution: L'Hospital rule was applied to $\frac{be^x -xe^b}{x-b}$. This implies that $g(x,b)$ is assumed to be continuous function but it is not at $x=b$. However, a finite value for $g(b,b)$ is obtained if $b \neq 0$.

Of course for either $x = 0$ or $b = 0$, $g(x,b)$ is infinite.

Hopefully, the above is a step in the right direction.