Re: Neumerical function



>>>>> "go" == go furuya <go_furuya@infoseek.jp> writes:

    go> Raymond toy write
    >> Footnotes:
    >> [1] This points out a small bug. maxima defines its own versions of
    >> asin, acos, atan, sinh, cosh, tanh, asinh, acosh, and atanh, which are
    >> all standard Lisp functions. This prevented me from computing
    >> elliptic functions for complex args. These should probably be
    >> modified. Plus I think the definitions for sinh, tanh, and perhaps
    >> others lose precision for x near 0.

    go> But I think your opnion wrong at sevral points.

Which points?  I am interested in your opinions.

    go> Especially at saying lose precision,you intend to say round errors,
    go> maxima treats fine this with bfloat.
    go> you shoud decide FPPREC length according to your goal precision.
    go> example

    go> (C1) FPPREC:100;
    go> (C2) sinh(1.0B-20),bfloat;
    go> (D2) 1.00000000000000000000000000000000000000001666666666666666666666666666666#

    go> 6666666666664063064574260973B-20

But consider this:

(C5) sinh(1.0d-4);

(D5) 				    1.0E-4

But sinh(x) = x + x^3/3! + ..., so the real answer should be closer to

(C6) 1.0d-4 + (1.0d-4^3)/6;

(D6) 			     1.000000001666667E-4

    go> this calculation is executed by according to source coad in
    go> trigi.lisp.
    go> (C3)exp(1.0B-20);
    go> (D3) 1.00000000000000000001000000000000000000005000000000000000000016666666666#

    go> 6666666667083333333333333333B0
    go> (C4)exp(-1.0B-20);

    go> (D4) 9.99999999999999999990000000000000000000049999999999999999999833333333333#

    go> 33333333375B-1
    go> (C5) D3-D4;
    go> (D5) 2.00000000000000000000000000000000000000003333333333333333333333333333333#

    go> 3333333333328126129148521946B-20
    go> (C6) %/2;
    go> (D6) 1.00000000000000000000000000000000000000001666666666666666666666666666666#

    go> 6666666666664063064574260973B-20

    go> We control round error with this mechanism and our head,I think.

Well, I don't think we're using our heads in this case.  

sinh(x) = 1/2*(exp(x) - exp(-x))
        = x + x^3/3! + x^5/5! + ...

So surely for small enough x, we should use our heads and use the
Taylor series (or something better) to evaluate sinh:

(C7) 1.0b-20 + (1.0b-20^3)/6;

(D7) 1.00000000000000000000000000000000000000001666666666666666666666666666666#

6666666666666666666666666667B-20


    go> But we shoud introduce "Interval Computation" in maxima ,like maple or
    go> mathematica.

Note that CMUCL has a interval arithmetic package that is used by the
compiler to derive the types of numbers.  It's correct for integers,
but for floating-point, it doesn't round up and down appropriately.

Ray