[Octave-bug-tracker] [bug #52074] Missing integration functions integral2, integral3

[David Bateman]

Follow-up Comment #12, bug #52074 (project octave):

Why mask at all ? You just need to invert the order of the integration in
dblquad and adapt the inner integrand to take into account the parametric y
boundaries, which is pretty easy. Basically as integral2 deals with the
argument checking, a couple of lines  imported from dblquad is enough for
this.

As for the 'tiled' version of the code, the Matlab version of this code uses
the TwoD function developed by Shampine and available at 

http://faculty.smu.edu/shampine/current.html

The license of Shampine is not clear though the code contains a clear
copyright line. So read the article but not the code!!!

However, the R developers have reimplemented the intaregal2 and integral3
functions based on Shampine's work in the "pracma" package available at

https://cran.r-project.org/web/packages/pracma/index.html

under a GPLv3 license. Converting the R code to Octave seems relatively
trivial, so see the attached function that I converted. It does the following

* Adds the parametric 'integrated' y boundaries
* Ports the R code to Octave
* Tries to auto detect the edge singularities in the 'tiled' method and
restarts the integration with singularity weakening transforms
* Adds the possibility to use unvectorized functions by wrapping them in array
fun

As an example of the performance of the new method 


>> format long
>> tic; integral2 (@(x, y) 1 ./ (x + y), 0, 1, 0, @(x) 1 -a x, "method",
"iterated"), toc
ans =  0.999999999999152
Elapsed time is 15.6569 seconds.
>> tic; integral2 (@(x, y) 1 ./ (x + y), 0, 1, 0, @(x) 1 - x, "method",
"tiled"), toc
ans =  0.999999962108304
Elapsed time is 0.172009 seconds.


This seems to me to be better version of integral2 and certainly more
complete. Though as I've done very little Octave coding for a number of years,
having someone read over the code and test it is probably a good idea. I also
took very little care with the Octave coding guidelines, so Rik go ahead and
rip my code apart ;-)

I also don't have a mercurial copy of Octave on my machine, so someone else
will have to work this up as a changeset.

D.

(file #41949)
    _______________________________________________________

Additional Item Attachment:

File name: integral2.m                    Size:19 KB


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?52074>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/