for loop results in NaN after 453 loops

I am trying to approximate an infinte sum expression. I dicovered that for whatever reason the sum becomes NaN after 453 loops.

Here is a sample of my code (just the portion I am having trouble with):

L = 1; H = 1; T0 = 373; T1 = 473;

theta_center = 0;

for p = 1:452 theta_center = theta_center + ((-1)^(p+1)+1)/p*sin(p*pi/2)*sinh(p*pi/L*H/2)/sinh(p*pi*H/L);

end

When I set the for loop to be any number less than 453, theta_center is a number, but when I set the for loop to any number greater 452, theta center becomes NaN.

Any help?

Answer 1

[Geoff]

Because the result of the first sinh term becomes Inf:

> p = 452;
> sinh(p*pi/L*H/2)

ans =

    1.1169e+308

> p = 453;
> sinh(p*pi/L*H/2)

ans =

    Inf

Note that the second sinh() term becomes Inf much sooner, but you only get NaN once you have divided Inf by Inf.

Do you want to treat p*pi as cyclic? Then use:

mod(p,2)*pi

EDIT: from comments below

I generated each iteration of your loop into a vector and then did a cumulative sum.

vals = [];
for p = 1:452
  vals(p) = ((-1)^(p+1)+1)/p*sin(p*pi/2)*sinh(p*pi/L*H/2)/sinh(p*pi*H/L);
end
vals = cumsum(vals);

% What is the last non-zero index?
idx = find( vals ~= 0, 1, 'last' );

That finds the real zeros, but it will be effectively zero much sooner:

% What is the last near-zero index?
idx = find( abs(vals) > 1e-8, 1, 'last' );

Without hard-coding any loop sizes, you can just use a while-loop with some terminating criteria: ie the most recent non-zero term being added to theta_center is sufficiently small.

[Phillip]

At what point point does matlab out put Inf? My TI-89 gives a value of 5.77375e+617 for the first sinh @ p=453. It also gives 0. for the entire term at p=453.

Also, i'm not sure what you mean by treating p*pi as cyclic.

[Geoff]

Your computer performs calculations using 64-bit double-precision numbers. The maximum representable value is approximately 1.79e+308. You might need to try expanding the sinh() series using various trig identities and see if you can keep your intermediate terms below the numerical limits of doubles.

Ignore the comment about p*pi being cyclic.

Anyway, with the parameters you gave, the solution becomes truly zero after 225 iterations, and effectively zero well before then. Why do you need to keep iterating past that point?

[Malcolm Lidierth]

Try

>> x=5.77375e+617

x =

Inf

See e.g. http://www.mathworks.com/company/newsletters/news_notes/pdf/Fall96Cleve.pdf

[Phillip]

Apparently I don't need to go past 225 iterations. I was unsure of how many iterations i needed, so I just tried 1000 to start with. And when I got NaN I became curious.

Btw, how did you determine that it becomes truely zero at 225 iterations?

FYI, I am doing a Finite Difference Approximation on a 2-D steady-steady heat transfer problem. After playing with the numbers I decided 10 iterations was plenty.

Thanks for all the help!

[Geoff]

Have edited my answer to explain. Would write it here, but the comments don't allow code-formatting.