Can't calculate e^(-x) for "high" x

Question

This is my code:

def fact(y):
    if y == 0:
        fact=1
        return fact
    else:
        fact=1
        for k in range (1, y+1):
            fact = fact * k
        return fact

def e_negative_x(x):
    n=0
    numerical_precesion=1
    numerical_precesion_ideal= 10**(-8)
    while numerical_precesion > numerical_precesion_ideal:
        sum=0
        for i in range (0, n+1):
            ind=((-x)**i)/fact(i)
            if i> 0 & n-i == 0:
                ind_2=((-x)**(i-1))/fact(i-1)
                numerical_precesion_1 = ind_2 - ind
                numerical_precesion = abs(numerical_precesion_1)
            sum = sum + ind
        if   numerical_precesion > numerical_precesion_ideal:
            n += 1
        elif sum < 0:
            n += 1
    return sum

I tried to use this for x=0.1;1;10;15 as a exercise and I got "correct" for my precision but when I tried x=30 it got stuck at -8e-5 (incorrect answer). I tried to increase my terms but it was still stuck at -8e-5 with a different decimal places and after a big increase in terms didn't change at all.

Edit: It was -8e-5 which is wrong because at infinity this series tend to 0. I made a infinite loop (for infinite terms) and made it print out how many terms it had and the sum for that terms. I had 600 + terms before i shut down. At 89term i got stuck at -8e5 and after 117 i got stuck at -8.553016433669241e-05 till 600 + term.

Doesn't seem it got stuck to me. I tried for 50 and it resulted as -7015.776232596799, whereas for 100 it resulted as -2.8756582514726483e+26. Are you sure these are not expected results? — mulaixi, Oct 20 '20 at 08:14
@mulaixi But these are obviously wrong values. The real function monotonically decreases towards 0 and is always positive. — Konrad Rudolph, Oct 20 '20 at 08:20
I didn't think about the function, yeah you're right. When I try the code the results fluctuate. So, these downvotes don't make sense. @KonradRudolph I didn't downvote btw. — mulaixi, Oct 20 '20 at 08:26
Unrelated, but you are shadowing the `sum` builtin and your own `fact` function inside of itself; also, you can use `math.factorial`. Also, your `elif` seems pointless: If the previous `if` is false, the loop will end and `n` is never used again. — tobias_k, Oct 20 '20 at 10:53
Also, `if i> 0 & n-i == 0:` does not do what you expect. This is not `(i> 0) and (n-i == 0)` but `i> (0 & n-i) == 0`, with `&` being bitwise-and, which boils down to just `i > 0` — tobias_k, Oct 20 '20 at 10:56
@tobias_k what u mean with shadowing my `sum` and `fact`? Actually I can't, it is asked to begin "new": no prior function. oh, ok i thought with putting elif i was putting another condition but i knew it was not working, how could i put it? `if i> 0 & n-i == 0:` It was doing fine, but i changed for `(if i> 0) & (n-i == 0):` anyway. — Miguel Godinho, Oct 20 '20 at 11:09
Also, you do not need two nested loops. Seems like you compute the entire taylor series up to n _again_ instead of just adding one more term to the last one. Does not fix the problem though. — tobias_k, Oct 20 '20 at 11:16

Lukas S · Accepted Answer · 2020-10-21T20:13:59.967

How to do it better:

Well I would say summing up positive and negative terms of comparable size is prone to error since it leads to summing small and big numbers which is known to be problematic. But you can avoid it since exp(-30)=1/exp(30).

1/math.fsum(((30)**n/math.factorial(n) for n in range(0,100+1)))

is

9.357622968840174e-14

i.e. positve and small like you would expect. And

1/math.fsum(((30)**n/math.factorial(n) for n in range(0,100+1)))-np.exp(-30)

gives

-1.262177448353619e-29

So we get bascially the same as numpy does.

Where exactly is the error:

Here is a table of different ways to calculate exp(-30) approximately with the taylor series summing up N terms:

Naive menas I just calculated the terms as floats and sumed them up.
In the fractions column I did the same but used python's fraction instead of float division.
Then I thought it would be nice to distinguish error in division/muliplication and in summing up. So I calculated each summand with fraction but turned them into floats before summing up.
Finally I calculated the values of exp(30) nativly and then 1/exp(30) which I am calling the 1/e^30_trick.

            native     fractions  float(fraction)  1/e^30_trick
0    1.000000e+00  1.000000e+00     1.000000e+00  1.000000e+00
10   1.212548e+08  1.212548e+08     1.212548e+08  4.187085e-09
20   8.529171e+10  8.529171e+10     8.529171e+10  2.652040e-12
30   3.848426e+11  3.848426e+11     3.848426e+11  1.706501e-13
40   6.333654e+10  6.333654e+10     6.333654e+10  9.670058e-14
50   8.782292e+08  8.782292e+08     8.782292e+08  9.360412e-14
60   1.685584e+06  1.685584e+06     1.685584e+06  9.357627e-14
70   6.225246e+02  6.225247e+02     6.225246e+02  9.357623e-14
80   5.588481e-02  5.595324e-02     5.588481e-02  9.357623e-14
90  -6.697346e-05  1.459487e-06    -6.697346e-05  9.357623e-14
100 -6.843293e-05  1.276217e-11    -6.843293e-05  9.357623e-14
110 -6.843294e-05  9.361706e-14    -6.843294e-05  9.357623e-14
120 -6.843294e-05  9.357623e-14    -6.843294e-05  9.357623e-14

I was quite happy with the result. Since it shows both the naiv and even the float(fractions) way are terrible so the error must be in summing up the negative and positive terms. Also even though both the 1/exp(30) trick and the factions do converge to the "correct" value of 9.357623e-14. The former converges much faster than the latter.

Table Code:

series = pd.Series(np.arange(0,151,10))

@np.vectorize
def naiv(N):
    return math.fsum(((-30)**n/math.factorial(n) for n in range(0,N+1)))

@np.vectorize
def fractions(N):
    return float(sum(Fraction((-30)**n,math.factorial(n)) for n in range(0,N+1)))

@np.vectorize
def float_fractions(N):
    return math.fsum(float(Fraction((-30)**n,math.factorial(n))) for n in range(0,N+1))

@np.vectorize
def one_over_30_trick(N):
    return 1/math.fsum(((30)**n/math.factorial(n) for n in range(0,N+1)))

pd.DataFrame({"native":naiv(series),"fractions":fractions(series),"float(fraction)":float_fractions(series),"1/e^30_trick":one_over_30_trick(series)},index=series)

Can't calculate e^(-x) for "high" x

1 Answers1

How to do it better:

Where exactly is the error: