The old round off problem?

[sam]

Hello all, I am taking a class in scientific programming at the local
college. My problem is that the following difference produces round off
errors as the value of x increases. For x >= 19 the diference goes to
zero.I understand the problem, but am curious as to whether their
exists a solution. I have tried various graphing programs,and they all
exihibit this problem.

thanks in advance Sam Schulenburg
f(x) = cosh^2(x) - sinh^2(x) = 1

from math import *
for x in range(20): print "x= %2d Sinh^2(x) = %20.3f f(x) = %2.10f
"%(x,pow(cosh(x),2),pow(cosh(x),2)- pow(sinh(x),2))
x= 0 Sinh^2(x) = 1.000 f(x) = 1.0000000000
x= 1 Sinh^2(x) = 2.381 f(x) = 1.0000000000
x= 2 Sinh^2(x) = 14.154 f(x) = 1.0000000000
x= 3 Sinh^2(x) = 101.358 f(x) = 1.0000000000
x= 4 Sinh^2(x) = 745.740 f(x) = 1.0000000000
x= 5 Sinh^2(x) = 5507.116 f(x) = 1.0000000000
x= 6 Sinh^2(x) = 40689.198 f(x) = 1.0000000000
x= 7 Sinh^2(x) = 300651.571 f(x) = 0.9999999999
x= 8 Sinh^2(x) = 2221528.130 f(x) = 1.0000000000
x= 9 Sinh^2(x) = 16414992.784 f(x) = 1.0000000037
x= 10 Sinh^2(x) = 121291299.352 f(x) = 1.0000000298
x= 11 Sinh^2(x) = 896228212.033 f(x) = 0.9999998808
x= 12 Sinh^2(x) = 6622280532.961 f(x) = 1.0000019073
x= 13 Sinh^2(x) = 48932402357.710 f(x) = 0.9999923706
x= 14 Sinh^2(x) = 361564266073.369 f(x) = 0.9999389648
x= 15 Sinh^2(x) = 2671618645381.616 f(x) = 1.0000000000
x= 16 Sinh^2(x) = 19740740045670.668 f(x) = 0.9921875000
x= 17 Sinh^2(x) = 145865435631864.219 f(x) = 1.0000000000
x= 18 Sinh^2(x) = 1077807886778799.250 f(x) = 1.0000000000
x= 19 Sinh^2(x) = 7963982939278438.000 f(x) = 0.0000000000

[Paul Rubin]

"sam" <sa******@pacbell.net> writes:

Hello all, I am taking a class in scientific programming at the local
college. My problem is that the following difference produces round off
errors as the value of x increases. For x >= 19 the diference goes to
zero.I understand the problem, but am curious as to whether their
exists a solution. [f(x) = cosh^2(x) - sinh^2(x) = 1]

The usual way is to notice that the difference goes to zero because
the lowest order terms of the Taylor series for cosh^2 and sinh^2 are
equal and cancel each other out. The solution is to write down the
series for (cosh^2(x) - sinh^2(x)) and add up a few non-cancelled terms.
All these series should converge very fast.

[David Treadwell]

On Mar 4, 2006, at 4:33 PM, Paul Rubin wrote:

"sam" <sa******@pacbell.net> writes:

Hello all, I am taking a class in scientific programming at the local
college. My problem is that the following difference produces
round off
errors as the value of x increases. For x >= 19 the diference goes to
zero.I understand the problem, but am curious as to whether their
exists a solution. [f(x) = cosh^2(x) - sinh^2(x) = 1]

The usual way is to notice that the difference goes to zero because
the lowest order terms of the Taylor series for cosh^2 and sinh^2 are
equal and cancel each other out. The solution is to write down the
series for (cosh^2(x) - sinh^2(x)) and add up a few non-cancelled
terms.
All these series should converge very fast.

Here's my analysis:

First, keep in mind that the range of values for a Taylor expansion
of sinh() and cosh() is limited. Then...

Write down the series for (cosh^2(x) - sinh^2(x)), as Paul suggested.

c2 = taylor(cosh(x)^2,x,8)

s2 = taylor(sinh(x)^2,x,8)

The first method gives the expected answer:

eval(c2-s2) ==> 1

Write down the series for cosh(x), square that; write the series for
sinh(x), square that and subtract.

Compare the two (I just did this in Eigenmath). Higher-order terms
remain for the second method.. I don't know how many terms my cosh()
and sinh() functions use, but if I evaluate (again, in Eigenmath)

The second method, which is more akin to what Python is doing (unless
there is an explicit 'sinh**2(x)' function) is different.

c2 = (taylor(cosh(x),x,8))^2, which gives terms up to and including
(x**8)**2 = x**16

s2 = (taylor(sinh(x),x,9))^2, which gives terms up to and including
x**18, then subtract the two, the result is

eval(c2 - s2) ==> 1 - 1/1814400 * x**10 - ... higher terms.

This begs the question, "What is the algorithm used by my Python?"

And, the ever-important question, "Does it matter in my final
result?" This is important, because I do not, in general, trust the
result when subtracting two large numbers.

:--David T.

[sam]

David I beg I beg
Can you answer the question?

Also thanks for the information on using the Taylor series.

Sam Schulenburg

[David Treadwell]

I wish I knew!

So I asked Google. Here's what I learned:

Most implementations are based on, or similar to the implementation
in the fdlibm package.

sinh(x) and cosh(x) are both based on exp(x). See http://
http://www.netlib.org/cgi-bin/netlib...dlibm/e_sinh.c

exp(x) is implemented by:

1. reducing x into the range |r| <= 0.5 * ln(2), such that x = k *
ln(2) + r
2. approximating exp(r) with a fifth-order polynomial,
3. re-scaling by multiplying by 2^k: exp(x) = 2^k * exp(r)

sinh(x) is mathematically ( exp(x) - exp(-x) )/2 but it's not
calculated directly that way.

My brain hurts at this point; it's late. I'll have another go at
hunting down the errors tomorrow. In the interim, does anybody out
there already know the answer?

:--David

On Mar 4, 2006, at 8:27 PM, sam wrote:

David I beg I beg
Can you answer the question?

Also thanks for the information on using the Taylor series.

Sam Schulenburg

--
http://mail.python.org/mailman/listinfo/python-list

[sam]

David Treadwell wrote:

exp(x) is implemented by:

1. reducing x into the range |r| <= 0.5 * ln(2), such that x = k *
ln(2) + r
2. approximating exp(r) with a fifth-order polynomial,
3. re-scaling by multiplying by 2^k: exp(x) = 2^k * exp(r)

sinh(x) is mathematically ( exp(x) - exp(-x) )/2 but it's not
calculated directly that way.

My brain hurts at this point; it's late. I'll have another go at
hunting down the errors tomorrow. In the interim, does anybody out
there already know the answer?

:--David

I tried the exp(x) equivalent of sinh(x) and cosh(x) with the same
results.
I think I will write my own or steal the Taylor(f(x),x,n) function in
both C++, and python.