bugGNU Octave - Bugs: bug #59403, integral2 and quad2d functions do...

 
 

bug #59403: integral2 and quad2d functions do not accurately document form of integrand function F

Submitter:  Luis Kao <kogiokka>
Submitted:  Tue 03 Nov 2020 06:15:14 AM UTC
   
 
Category:  Documentation Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Luis Kao Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 14 Feb 2024 03:49:04 AM UTC, comment #14: 

All seems good here.  Marking as Fixed and closing report.

Rik <rik5>
Group administrator
Wed 31 Jan 2024 05:30:21 AM UTC, comment #13: 

pushed a docstring update for quad2d, integral2, integral3 as
https://hg.savannah.gnu.org/hgweb/octave/rev/2da9214df7a4

will leave report open for now in case there's any desire to look for better ways to catch input errors to these functions. if not, we can probably close this as fixed and consider the doc warning sufficient.

Nicholas Jankowski <nrjank>
Group Member
Wed 31 Jan 2024 03:27:11 AM UTC, comment #12: 

the test cases for integrands with non-vectorized inputs:

1d.  The help for all of the 1D integral functions (except quadv) clearly state something like "must be vectorized and return a vector of output values when given a vector of input values". quadv or 'arrayvalued' options force the integrator to provide scalar inputs and avoid matrix math problems.


>>f = @(x) sinc(x) * sinc(x);

>> quadv(f, -1, 1)  ## designed for vector/matrix valued fns
ans = 0.9028

>> integral(f, -1, 1)  ## uses quadcc
error: operator *: nonconformant arguments (op1 is 33x1, op2 is 33x1)
error: called from
    integral at line 146 column 11

>> quadgk(f,-1,1)
error: operator *: nonconformant arguments (op1 is 150x1, op2 is 150x1)
error: called from
    @<anonymous> at line 1 column 18
    quadgk>@<anonymous> at line 355 column 28
    quadgk>__quadgk_eval__ at line 558 column 5
    quadgk at line 369 column 22

>> integral(f, -1, 1, 'arrayvalued',true)
ans = 0.9028

>> quadgk(f, -1, 1, 'arrayvalued',true)
ans = 0.9028


2D:
integral2, dblquad, and quad2d state: "F must be of the form z = f(x,y) where X is a vector and Y is a scalar.  It should return a vector of the same length and orientation as X."  As per last comment, this is actually only valid for dblquad (and integral2 when it mimics dblquad with the 'iterated' method). 'vectorized' option for integral2 & dblquad wrap f with arrayfun to force scalar operations in the integrand.


>> f = @(x, y) sinc(x) * sinc(y);

>> integral2(f, -1, 1, -1, 1)
ans = 12.328

>>dblquad(f, -1, 1, -1, 1)
ans = 1.3900

>> quad2d(f, -1, 1, -1, 1)
ans = 12.328

>> integral2(f, -1, 1, -1, 1, 'vectorized', true)
ans = 12.328

>> integral2(f, -1, 1, -1, 1, 'vectorized', false)
ans = 1.3900

>> integral2(f, -1, 1, -1, 1, 'method', 'iterated') #avoids quad2d
ans = 1.3900

>> quad2d(f, -1, 1, -1, 1, 'vectorized', true)
ans = 12.328

>> quad2d(f, -1, 1, -1, 1, 'vectorized', false)
ans = 1.3900


matlab produces an input error for integral2 and quad2d for having non-vectorized operations, but dblquad produces the same output as octave.

3D:

>> f = @(x, y, z) sinc(x) * sinc(y) * sinc(z);

>> integral3(f, -1, 1, -1, 1, -1, 1) ##uses quad2d
ans = 14.535

>> triplequad(f, -1, 1, -1, 1, -1, 1)
ans = 1.6388

>> integral3(f, -1, 1, -1, 1, -1, 1, 'method', 'iterated')
ans = 1.6388

>> integral3(f, -1, 1, -1, 1, -1, 1, 'vectorized', true)
ans = 14.535

>> integral3(f, -1, 1, -1, 1, -1, 1, 'vectorized', false)
ans = 1.6388


(where matlab again errors for integral3 and produces the same answer for triplequad)

Nicholas Jankowski <nrjank>
Group Member
Wed 31 Jan 2024 03:17:53 AM UTC, comment #11: 

I didn't really give this much attention a few years ago.  (trying to graduate, etc.)

reading through things to see where a doc change might be made, I'm think that quad2d may actually be doing something wrong when it's passed a not-properly-vectorized function (as opposed to it just doing a proper integral on the non-elementwise function it's given). for this f, only quad2d and the integral2 and 3 functions calling quad2d are generating different numerical outputs without erroring. generally, it's important that integration query points be independent (otherwise the summation loses meaning), and sinc(x)sinc(y) still does that when y is a scalar. integral2 and 3 select between iterated (using quadcc), and quad2d. Iterated produces a correct result.  quad2d doesn't.

The comments below basically capture the situation. When I made the 1st stab at integral2, I retained the input requirements from dblquad (vector x, scalar y) as it was really just a wrapper.  then we (well david) improved it to actually do a decent job, creating quad2d, too, but we mistakenly kept those same input requirements for f in the doc.  As rik saw, stepping through quad2d, at the end of the tensorproduct subfunction it calls f(x,y) after expanding both x and y to  2D matrices like they're doing a manual 'meshgrid'. That appears in fact to be a feature of the Shampine algorithm our quad2d (and matlab's according to their help doc) is based on. This sort of requires elementwise operation so it can do vectorized expansion for speed. if the elements aren't independent, it loses meaning. iterated methods sometimes get away with odd functions, but by permitting non-elementwise operations, we're really violating the intent of the Shampine method. most of the integrators have options to turn off vectorized solving specifically to navigate around this problem.

right now 1d, 2d, and 3d integrators wind up failing with obscure errors for malformed f's, catching them would be better.  And apparently quad2d may have been silently producing erroneous results since octave 4. from jwe's comment #2, I'm guessing we don't have the ability to throw an upfront error? that would be ideal. the only thing I could think of is something crude like assert(vectorize(f),f).  but even vectorize says not to use vectorize. and that won't catch much anyway.

If the best we can do is a doc change, i think it needs to state tighter input requirements, specifically as elementwise.  "inputs need to have x form or numerical results can be unpredictable. if you can't then use 'vectorized=false'" that way at least people are warned, with a workaround, and silent failures like comment #0 can be waved away with 'we told you so'.  I'll generate some text and see how it reads.

Nicholas Jankowski <nrjank>
Group Member
Sat 14 Nov 2020 03:00:45 PM UTC, comment #10: 

comment #9:

> We can't base the implementation of Octave on any information obtained from looking at parts of Matlab that might be available.  If you value the existence of Octave, please do not peek at and share information about the internals of Matlab here.


I'm sorry. It won't happen again.

Luis Kao <kogiokka>
Sat 14 Nov 2020 02:12:22 PM UTC, comment #9: 

We can't base the implementation of Octave on any information obtained from looking at parts of Matlab that might be available.  If you value the existence of Octave, please do not peek at and share information about the internals of Matlab here.

John W. Eaton <jwe>
Group administrator
Sat 14 Nov 2020 09:05:21 AM UTC, comment #8: 

comment #7:

> ... The comment states that it would make it harder for the user "generating an error by using matrix functions instead of elementwise functions". ...


EDIT:
otherwise it would make it harder for the user

Luis Kao <kogiokka>
Sat 14 Nov 2020 09:00:34 AM UTC, comment #7: 

So I took a peek at quad2d in MATLAB. It does check if "FUN is properly vectorized". The comment states that it would make it harder for the user "generating an error by using matrix functions instead of elementwise functions". I don't know if there is a reason using matrix function within integration, but MATLAB seems to suggest users should not do that.

Also there is another thing they explicitly check about size mismatch:


quad2d(@(x,y) 1, 0, 1, 0, 1)


Octave returns value 1, while MATLAB throws an error: "Integrand output size does not match the input size."

Luis Kao <kogiokka>
Thu 12 Nov 2020 02:03:10 PM UTC, comment #6: 

Matlab uses a G3,K7 Guass-Konrad rule and Octave a G7,K15 Guass-Konrad rule in quad2d. So the only difference is the matlab should call the integrand with a 7x7 matrix of scalar arguments and Octave with a 15x15 matrix of scalar arguments. So matlab is detecting the issue somehow. Perhaps they do something like


try
   f(ones(2,1),ones(2,1)
catch
   error('ERROR HERE')
endcatch


To detect matrix cases before the rest of quad2d. But if they are then calling the integrand better not have any side effects

D.


David Bateman <dbateman>
Group Member
Thu 05 Nov 2020 05:01:20 PM UTC, comment #5: 

This helps make my argument.  The documentation for the form of the integrand function is misleading.  It used to be the case that F had to take the very particular form mentioned in integral2/quad2d because of the way 2-D integration was done.  This is no longer the case, but you can still see a bit of the history of the function in that there is a "method" property that can take the value "iterated" or "tiled".  For "iterated", F used to have that specific form.

Rik <rik5>
Group administrator
Thu 05 Nov 2020 04:25:08 PM UTC, comment #4: 

After seeing the comments and checking the documentation again, I would like to provide some extra arguments.




If the user's intention is to do matrix multiplication, they will expect a matrix output. For example


function z = f(x,y)
  z = sinc(x) * sinc(y);
endfunction

x = [1 2; 3 4];
y = [5 6; 7 8];

f(x, y)
    =>   2.0836e-49  -3.7300e-50
         2.0836e-49  -3.7300e-50


And since that's the form of the original function, I suspect that they would want a matrix output from integral2/quad2d as well.


Another thing is that the documentation says

>     F is a function handle, inline function, or string containing the
>     name of the function to evaluate.  The function F must be of the
>     form z = f(x,y) where X is a vector and Y is a scalar.  It should
>     return a vector of the same length and orientation as X.


I don't know the implementation reason why Y should be a scalar, but it seems like current behaviors of integral2/quad2d disagree with the documentation.


function z = f(x,y)
  disp("Y isscalar:");
  isscalar(y)
  z = sinc(x) * sinc(y);
endfunction

integral2(@f, -1, 1, -1, 1);


The result from the script above is

Y isscalar:
ans = 0
Y isscalar:
ans = 0
Y isscalar:
ans = 0
Y isscalar:
ans = 0


Luis Kao <kogiokka>
Thu 05 Nov 2020 04:20:31 PM UTC, comment #3: 

This is what I suspected.  It was more or less luck, an artifact of how Matlab implements integral2, that led to the warning.  In Octave, I think the underlying implementation is trying for better performance so instead of calling the integrand function 15x15 = 225 times and incurring the overhead of a function call it calls the integrand once with a meshgrid 15x15 input.

Note, that Matlab's own documentation says the integrand function must use element-wise operators.


Integrand, specified as a function handle, defines the function to be integrated over the planar region xmin ≤ x ≤ xmax and ymin(x) ≤ y ≤ ymax(x). The function fun must accept two arrays of the same size and return an array of corresponding values. It must perform element-wise operations.


My leaning is to treat this bug as a documentation issue and more specifically mention that the function must return an array and use element-wise operators.

Rik <rik5>
Group administrator
Tue 03 Nov 2020 05:52:41 PM UTC, comment #2: 

I changed the example to be


function r = fxy (x, y)
  size (x), size (y),
  r = sinc (x) * sinc (y);
end
integral2 (@fxy, -1, 1, -1, 1)


and I see


ans =

   15   15

ans =

   15   15

ans =

   15   15

ans =

   15   15

ans =

   15   15

ans =

   15   15

ans =

   15   15

ans =

   15   15

ans = 12.328


So, unfortunately, the integration function has selected X, Y values that do conform for matrix multiplication.  So I don't see how we would detect an error here.

John W. Eaton <jwe>
Group administrator
Tue 03 Nov 2020 05:31:37 PM UTC, comment #1: 

Confirmed.  Changing release to "dev" as any fix is likely to occur on the development branch so there is enough time to thoroughly test a solution before release.

The solution may be more documentation.  Both functions are valid and therefore a user might have wanted to integrate the first one, even if it does seem like an odd choice.

Another possibility, if the condition can be detected easily, would be to emit a warning rather than an error so the user would have the chance to check whether the integrand function was really the one they intended.

Rik <rik5>
Group administrator
Tue 03 Nov 2020 06:15:14 AM UTC, original submission:  

Octave silently produces an "incorrect" result when the input function F does not use element-by-element multiplication. This behavior is different from MATLAB, which throws the error message:

```
Error using  *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To perform elementwise multiplication, use '.*'.
```

Example Script


f = @(x, y) sinc(x) * sinc(y);
integral2(f, -1, 1, -1, 1)
    => 12.32805347760951

f = @(x, y) sinc(x) .* sinc(y);
integral2(f, -1, 1, -1, 1)
    => 1.389993237875656

Luis Kao <kogiokka>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by dbateman (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by kogiokka (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2024-02-14 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2024-01-31 nrjank StatusConfirmed Ready For Test
    2020-11-05 rik5 CategoryOctave Function Documentation
        Item GroupMissed Error or Warning Inaccurate Result
        Operating SystemGNU/Linux Any
        Summaryintegral2 and quad2d functions do not catch incorrect dimensions for matrix multiplication integral2 and quad2d functions do not accurately document form of integrand function F
    2020-11-03 rik5 StatusNone Confirmed
        Release5.2.0 dev

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code