octave-bug-tracker
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Octave-bug-tracker] [bug #62321] Wrong result for Airy funtion


From: anonymous
Subject: [Octave-bug-tracker] [bug #62321] Wrong result for Airy funtion
Date: Mon, 18 Apr 2022 03:39:05 -0400 (EDT)

Follow-up Comment #14, bug #62321 (project octave):

Your patch in [comment #11] looks good to me.

Additionally, there is the "duplication" with the Airy function Bi, the
scaling is implemented in the zbiry
<https://hg.savannah.gnu.org/hgweb/octave/file/30f7f409861a/liboctave/external/amos/zbiry.f#l24>/cbiry
<https://hg.savannah.gnu.org/hgweb/octave/file/30f7f409861a/liboctave/external/amos/cbiry.f#l23>
FORTRAN code, but it is done in the ocatve code (l1388
<https://hg.savannah.gnu.org/hgweb/octave/file/30f7f409861a/liboctave/numeric/lo-specfun.cc#l1388>
and l1462
<https://hg.savannah.gnu.org/hgweb/octave/file/30f7f409861a/liboctave/numeric/lo-specfun.cc#l1462>).

As for some test, with your finding, we can reduce the bug to the simple
example:

## Just the first element interests us here, but we need to put it in a
complex vector,
## so it remembers that its `-1-0i` and not just `-1`.
z = [(-1-1i) + 1i, 1+1i];

## The error should be 0 for Ai
(airy(0, z) - airy(0, conj(z)))(1)

## and for Ai'.
(airy(1, z) - airy(1, conj(z)))(1)

and I have check that this code does not give machine precision on my `6.2.0`
octave version.

You could also test the relations:

z = -1 + 1i;

## Those should be machine precision
airy(0, z, true) - airy(0, z, false) * exp(2/3 * z * sqrt(z))
airy(1, z, true) - airy(1, z, false) * exp(2/3 * z * sqrt(z))
airy(2, z, true) - airy(2, z, false) * exp(-abs(real(2/3 * z * sqrt(z))))
airy(3, z, true) - airy(3, z, false) * exp(-abs(real(2/3 * z * sqrt(z))))




    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?62321>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/




reply via email to

[Prev in Thread] Current Thread [Next in Thread]