[Octave-bug-tracker] [bug #56665] Incorrect results with ppval and interp1, with complex input

[Chad Oldfield]

URL:
  <https://savannah.gnu.org/bugs/?56665>

                 Summary: Incorrect results with ppval and interp1, with
complex input
                 Project: GNU Octave
            Submitted by: chadoldfield
            Submitted on: Wed 24 Jul 2019 08:50:41 PM UTC
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Incorrect Result
                  Status: None
             Assigned to: None
         Originator Name: Chad Oldfield
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: 5.1.0
        Operating System: Microsoft Windows

    _______________________________________________________

Details:

The following is a copy of the command window, obtained right after freshly
starting Octave 5.1.0:

>> interp1([1 2], [2 4], 1.5 + 0.01i)
ans =  3.000000 - 0.020000i
>> pp = interp1([1 2], [2 4], 'pp');
>> ppval(pp, 1.5 + 0.01i)
ans =  3.000000 - 0.020000i
>> interp1([-2 2], [-4 4], 1.5 + 0.01i)
ans =   NA

The expected output for all of these has the opposite sign for the imaginary
part:
ans =  3.000000 + 0.020000i

This can be correctly obtained, for example, using:
>> p = polyfit([1 2], [2 4], 1);
>> polyval(p, 1.5 + 0.01i)

The background behind this example is that it is convenient to use the complex
step method to numerically calculate the derivative of a complicated function.
 Such a complicated function may include interpolation (among various other
steps and algorithms).  The method is described in the link below.  Much of
the utility of this is that it avoids some of the numerical errors involved
with finite differences.

http://mdolab.engin.umich.edu/content/guide-complex-step-derivative-approximation-0

Essentially, the derivative of the real-valued function f can be approximated
by:

dfdx = imag(f(x + 1i * step)) / step

where step is a small value (eg. 1e-20).  The expected complex part of the
results of all of the tests above is therefore 0.01 times the slope of the
line, or 0.02.

In the examples below, the first two (interp1 and ppval) may be due to the use
of ppval by interp1, and the use of the transpose operator at line 71 of
ppval.m:

  dx = (xi - x(idx))(:)';

It is possible that applying the array transpose would help:

  dx = (xi - x(idx))(:).';

In the third example, it is possible that the NA result is due to the use of
the complex magnitude in the comparison operators in line 380 of interp1.m:

    outliers = (xi < minx) | ! (xi <= maxx);  # this even catches NaNs

When using the complex step method, the real-valued function needs to be
written such that the same path is followed with a real vs. a complex input. 
This could be achieved, for example, with:

    outliers = (real(xi) < minx) | ! (real(xi) <= maxx);

Thanks for your help!




    _______________________________________________________

Reply to this item at:

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

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