4

Does anyone know how to make the following Matlab code approximate the exponential function more accurately when dealing with large and negative real numbers?

For example when x = 1, the code works well, when x = -100, it returns an answer of 8.7364e+31 when it should be closer to 3.7201e-44.

The code is as follows:

s=1
a=1;
y=1;
for k=1:40
    a=a/k;
    y=y*x;
    s=s+a*y;
end
s

Any assistance is appreciated, cheers.

EDIT:

Ok so the question is as follows:

Which mathematical function does this code approximate? (I say the exponential function.) Does it work when x = 1? (Yes.) Unfortunately, using this when x = -100 produces the answer s = 8.7364e+31. Your colleague believes that there is a silly bug in the program, and asks for your assistance. Explain the behaviour carefully and give a simple fix which produces a better result. [You must suggest a modification to the above code, or it's use. You must also check your simple fix works.]

So I somewhat understand that the problem surrounds large numbers when there is 16 (or more) orders of magnitude between terms, precision is lost, but the solution eludes me.

Thanks

EDIT:

So in the end I went with this:

s = 1;
x = -100;
a = 1;
y = 1;
x1 = 1;

for k=1:40
    x1 = x/10;
    a = a/k;
    y = y*x1;
    s = s + a*y;
end

s = s^10;
s

Not sure if it's completely correct but it returns some good approximations.

exp(-100) = 3.720075976020836e-044
s = 3.722053303838800e-044

After further analysis (and unfortunately submitting the assignment), I realised increasing the number of iterations, and thus increasing terms, further improves efficiency. In fact the following was even more efficient:

s = 1;
x = -100;
a = 1;
y = 1;
x1 = 1;

for k=1:200
    x1 = x/200;
    a = a/k;
    y = y*x1;
    s = s + a*y;
end

s = s^200;
s

Which gives:

exp(-100) = 3.720075976020836e-044
s = 3.720075976020701e-044

EJG89
  • 1,189
  • 7
  • 17
Billy Ray Valentine
  • 401
  • 5
  • 17

2 Answers2

7

As John points out in a comment, you have an error inside the loop. The y = y*k line does not do what you need. Look more carefully at the terms in the series for exp(x).

Anyway, I assume this is why you have been given this homework assignment, to learn that series like this don't converge very well for large values. Instead, you should consider how to do range reduction.

For example, can you use the identity

exp(x+y) = exp(x)*exp(y)

to your advantage? Suppose you store the value of exp(1) = 2.7182818284590452353...

Now, if I were to ask you to compute the value of exp(1.3), how would you use the above information?

exp(1.3) = exp(1)*exp(0.3)

But we KNOW the value of exp(1) already. In fact, with a little thought, this will let you reduce the range for an exponential down to needing the series to converge rapidly only for abs(x) <= 0.5.

Edit: There is a second way one can do range reduction using a variation of the same identity.

exp(x) = exp(x/2)*exp(x/2) = exp(x/2)^2

Thus, suppose you wish to compute the exponential of large number, perhaps 12.8. Getting this to converge acceptably fast will take many terms in the simple series, and there will be a great deal of subtractive cancellation happening, so you won't get good accuracy anyway. However, if we recognize that

12.8 = 2*6.4 = 2*2*3.2 = ... = 16*0.8

then IF you could efficiently compute the exponential of 0.8, then the desired value is easy to recover, perhaps by repeated squaring.

exp(12.8)
ans =
          362217.449611248

a = exp(0.8)
a =
          2.22554092849247
a = a*a;
a = a*a;
a = a*a;
a = a*a
          362217.449611249

exp(0.8)^16
ans =
          362217.449611249

Note that WHENEVER you do range reduction using methods like this, while you may incur numerical problems due to the additional computations necessary, you will usually come out way ahead due to the greatly enhanced convergence of your series.

  • Thanks Woodchips, not sure how I would apply this to the algorithm, especially if x = -100. – Billy Ray Valentine Apr 23 '12 at 15:51
  • exp(hpf(-100)) ans = 3.720075976020835962959695803863118337358892292376781967120613876663290475895815718157118778642281497e-44 –  Apr 23 '12 at 16:19
5

Why do you think that's the wrong answer? Look at the last term of that sequence, and it's size, and tell me why you expect you should have an answer that's close to 0.

My original answer stated that roundoff error was the problem. That will be a problem with this basic approach, but why do you think 40 is enough terms for the appropriate mathematical ( as opposed to computer floating point arithmetic) answer.

100^40 / 40! ~= 10^31.

Woodchip has the right idea with range reduction. That's the typical approach people use to implement these kinds of functions very quickly. Once you get that all figured out, you deal with roundoff errors of alternating sequences, by summing adjacent terms within the loop, and stepping with k = 1 : 2 : 40 (for instance). That doesn't work here until you use woodchips's idea because for x = -100, the summands grow for a very long time. You need |x| < 1 to guarantee intermediate terms are shrinking, and thus a rewrite will work.

John
  • 5,735
  • 3
  • 46
  • 62
  • 1
    Well, rounding is one way of looking at it, but really not. The pertinent phrase is "massive subtractive cancellation". There are VERY large negative value, and VERY large positive values as terms in the series. The two will in theory, if you have a large number of terms, eventually end up canceling out. In practice, this does not happen, because the limits of floating point arithmetic cause you to lose precision. This is why the process is known as massive subtractive cancellation. –  Apr 23 '12 at 13:53
  • I was rewriting as you commented. Funny. – John Apr 23 '12 at 13:55
  • Yes, good timing. I'll leave my comment in there, only because the massive subtractive cancellation concept is an important one. –  Apr 23 '12 at 14:01
  • The problem is that there is meant to be a simple fix for the issue. Unless of course the code HAS to be completely rewritten. I understand that when you exceed 16 digits of accuracy, you lose precision, so with large negative values of x (with this algorithm) or larger positive values of x the difference between terms is to great to compute accurately. The question is, can this be solved with a simple fix? Or should one reduce the terms and accept that it will only work by adding boundary conditions? – Billy Ray Valentine Apr 23 '12 at 14:21
  • You'll have to give us the exact statement of your homework problem in order to get a good answer, I don't know what "simple fix" you are referring to. Woodchips domain reduction comment is the simplest fix I know of. – John Apr 23 '12 at 15:13
  • 1
    @BillyRayValentine - Why do you think range reduction is not simple? Either of the schemes I show can be accomplished in only a few lines of extra code. Admittedly, there are clever variations that can add several dozen lines to gain extreme convergence when that is important, but the simple range reductions I describe are short to implement. TRY IT! You will learn much in doing so. –  Apr 23 '12 at 16:01
  • @woodchips - I'm not saying it's not simple, but it looks like you're using the exponential function explicitly. Thinking about it, are you saying reduce the input x by some factor, say 10 and then compute the value of -10 as the new x then manipulate that value somehow? – Billy Ray Valentine Apr 23 '12 at 16:19
  • Ok so just tried reducing by 10 and then raising the answer to the power of 10... and it works! Thanks man. It's all about these identities! – Billy Ray Valentine Apr 23 '12 at 16:25
  • 1
    It really IS about those identities. Test your scheme out for various numbers, both positive and negative. If you have a problem, then ask again. But really, just working with those identities will give what you need. –  Apr 23 '12 at 17:29