Why my code doesn't work?

I m trying to make a program that calculates cosh(x) using an approximation. The error between the calculated value and the real value must be less than e-10. If i try to run the code for small numbers or very large the results are correct, but for x=50 for example, the variable termSum becomes NaN and everything is ruined. What am i doing wrong?

clc;
clear;

x = input('Give x to calculate cosh(x): ');
disp(' ')

y1 = (exp(x)+exp(-x))/2;

terms = zeros(1, 1500);
err = 1;
i=0;

while err>10^(-10)

      terms(i+1) = (abs(x).^i)/factorial(i);
      termSum = sum(terms);
      ex = termSum^sign(x);
      y2 = (ex + ex^(-1))/2;

      err=abs(y2-y1);

      i=i+1;

end;

disp(['cosh(', num2str(x),')=', num2str(y2)])

Answer 1

[Guillaume]

I've not tested your code but to me it does not look that it's going to produce the right result for any value.

For a start, I don't understand why you use a series expansion of the exponential rather than a straightforward expansion of cosh directly.

edit: removed my own nonsense

You're getting NaN for x = 50 because at the 182th iteration, 50^182 is Inf, and as explained on the factorial page, the factorial of any number above 170 is also Inf, so for i = 182, factorial is Inf. Inf/Inf is Nan.

There's no point predeclaring 1500 terms since as stated you can't get a factorial above 170 terms. You may as well put a hard stop to your calculation when i reach 170.

[liakoyras]

Well, thanks for answering. I have to do the approximation on e^x, it is part of the question.

As for the two lines, they are indeed calculating cosh, for the values of x the code works it returns the same value as (exp(x)+exp(-x))/2. And that 1500 was out of frustration, i didn't know where my mistake was :p.

Despite all that, i found your answer very helpful, thanks a lot.

Do you have any other ways of approximating e^x to recommend (besides the series expansion)?

[Guillaume]

Gah! I'm the one talking nonsense. You are indeed calculating cosh correctly. The problem is that you're never going to reach your tolerance for such big cosh values.

>> cosh(50)
ans =  
   2.5924e+21

>> eps(cosh(50))  
ans =  
      524288

This means that the nearest doubles to cosh(50) are cosh(50) &#177 52488. That is the minimum precision that you could achieve. Nowhere near 1e-10. In fact, you won't be able to reach your precision for any cosh value above or equal to 2^19, since that where eps crosses above 1e-10.

Since you can never achieve your precision, your iteration will keep on going forever. As I've already explained after step 170, it will produce Inf because of the factorial and from step 182, NaN since the numerator is also Inf.

[liakoyras]

So you mean that there is no way to achieve that precision on matlab?

[Guillaume]

It's not just matlab, it's inherent to floating point numbers.

You may be able to achieve this precision with the symbolic toolbox but I suspect that it would slow your algorithm to a crawl. Storing a number of 1e21 magnitude with a precision of 1e-10 is kind of meaningless.

You really need to adjust your precision to the magnitude of the result.

[liakoyras]

Thanks a lot, unfortunately I have to deal with this precision.

Answer 2

[Kelly Kearney]

I haven't given any thought to the algorithm you show, though a quick test does seem to indicate it converges to cosh for small numbers. But as Guillaume points out, both the numerator and denominator in your terms equation get pretty big pretty quickly:

nterm = 200;
ii = 1:nterm;

x = 50
tnum = abs(x).^ii;
tden = factorial(ii);

terms = tnum./tden;

You can examine the values here:

[tnum' tden' terms']

The denominator becoming infinite (after 172 terms, due to limitations in factorial) isn't actually a problem, although at that point the extra term adds 0 to the sum ( x/Inf = 0 ) and the convergence stops. If the numerator becomes infinite first (e.g. x = 1000), you'll start getting Inf in your terms ( Inf/x = Inf ), and if both become infinite, you'll get NaNs ( Inf/Inf = NaN ).