Vtrap derived

[JimH]

I have playing around with enough exponential functions lately that I started to wonder where the "vtrap" equation I was using came from.

The general form for this equation is:

vtrap = x/(exp(x/y) - 1)

however as can be easily seen if x/y is 0, or close to zero, then the denominator of the given equation is really small or zero, leading to and infinite or very large output. To get around this when x/y is really small the output is given as:

vtrap = y*(1 - x/y/2) or as I prefer y - x/2

hence the following code:

Code: Select all

FUNCTION vtrap(x,y) {  :Traps for 0 in denominator of rate eqns.
        if (fabs(x/y) < 1e-6) {
                vtrap = y*(1 - x/y/2)
        }else{
                vtrap = x/(exp(x/y) - 1)
        }
}

but how to go from one formula to the other?

To start exp(x/y) can be approximated by (where u = x/y) 1 + u + u^2/2 + more terms which make the approximation more accurate which I will not use (second order approximation). Thus from the first equation we have:

x/(u+u^2/2) => the 1's canceled

Going back to x/y the equation can be simplified fairly easily to:

y^2/(y+x/2)

however this still contains a denominator, which I guess would be nice to get rid of (??) or perhaps the squared term for y isn't stable (??) (if we assume that y may be very large to cause x/y to be small)
Regardless, the following simplification can be made:

y^2 * (y-x/2)
-------------------------
(y+x/2) * (y-x/2)

This leads to:

y^3 - y^2*x/2
-----------------
y^2 - x^2/4

If x^2/4 is assumed to be small, or I guess small compared to y^2, the denominator can be approximated as y^2 (this is the second and final approximation).

This simplifies the original equation down to y - x/2

[ted]

A shorter derivation starts as yours did, but branches off at this point:

x/(u+u^2/2) => the 1's canceled

Rewrite this as
(x/u) / (1 + u/2) = y / (1 + u/2)
and then take advantage of the fact that, for |epsilon| << 1,
1 / 1 + epsilon
is very closely approximated by
1 - epsilon

Thus
y / (1 + u/2) ~ y (1 - u/2) = y (1 - x/y/2)