# help with double quadrature

 Classic List Threaded
6 messages
Reply | Threaded
Open this post in threaded view
|

## help with double quadrature

 I have to evaluate an integral over a triangular area:int_0^1 int_0^y [f(x,y)] dx,dyInitially I though I could use dblquad, but is seems the integration limits have to be numbers, not variables.  I have backed of to using quadgk to do the inner integral, then summing over all y.  Does anyone know a better way? _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

## Re: help with double quadrature

 On Mon, Sep 11, 2017 at 8:17 PM, Clinton Winant <[hidden email]> wrote: > I have to evaluate an integral over a triangular area: > int_0^1 int_0^y [f(x,y)] dx,dy > > Initially I though I could use dblquad, but is seems the integration limits > have to be numbers, not variables.  I have backed of to using quadgk to do > the inner integral, then summing over all y.  Does anyone know a better way? Have you consider using Monte Carlo integration? You can use rejection sampling (you'll get about 50% in the unit square) or the hit-and-run sampler[1] The rejection samplers goes something like (assuming your triangle is defined for (x,y) in [0,1]x[0,1]) triag_boundary =@(x) 1 - x # your triangle N = 1e6; xy = rand(N,2); xy (xy(:,2) > triag_boundary(xy(:,1)),:) = []; I = sum (f(xy(:,1), xy(:,2))) / N; # this is your integral [1]: https://pdfs.semanticscholar.org/4dfe/8bde57bfe0d27bf48ed4d8e32f4938699392.pdf_______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

## Re: help with double quadrature

 On Mon, Sep 11, 2017 at 4:04 PM, Juan Pablo Carbajal wrote:On Mon, Sep 11, 2017 at 8:17 PM, Clinton Winant <[hidden email]> wrote: > I have to evaluate an integral over a triangular area: > int_0^1 int_0^y [f(x,y)] dx,dy > > Initially I though I could use dblquad, but is seems the integration limits > have to be numbers, not variables.  I have backed of to using quadgk to do > the inner integral, then summing over all y.  Does anyone know a better way? Have you consider using Monte Carlo integration?Never saw that before.  That's interesting.  How much better is in than, say, a 3D equivalent of 'trapz'? (is there a name for that? patch summation over a 2D mesh of points?In any case, yes, dblquad only works over a rectangular region.  Mathworks implemented variable limits into both quad2d and integral2, but Octave has not yet implemented compatible functions for those.Octave documentation does provide a possible double integration solution using collocation, I don't know if this is as easy to adapt to a non-rectangular region as suggested by Juan Pablo.https://www.gnu.org/software/octave/doc/v4.2.1/Functions-of-Multiple-Variables.html _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

## Re: help with double quadrature

 On Mon, Sep 11, 2017 at 3:31 PM, Nicholas Jankowski wrote:On Mon, Sep 11, 2017 at 4:04 PM, Juan Pablo Carbajal wrote:On Mon, Sep 11, 2017 at 8:17 PM, Clinton Winant <[hidden email]> wrote: > I have to evaluate an integral over a triangular area: > int_0^1 int_0^y [f(x,y)] dx,dy > > Initially I though I could use dblquad, but is seems the integration limits > have to be numbers, not variables.  I have backed of to using quadgk to do > the inner integral, then summing over all y.  Does anyone know a better way? Have you consider using Monte Carlo integration?Never saw that before.  That's interesting.  How much better is in than, say, a 3D equivalent of 'trapz'? (is there a name for that? patch summation over a 2D mesh of points?https://en.wikipedia.org/wiki/Monte_Carlo_integration​It is cute, but I found that it is not very practical if you want some decent precision.​​(But I have not tried some advanced algorithms there.)​Dmitri.​ _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

## Re: help with double quadrature

 Matlab's documentation suggests using dblquad with the function modified to be 0 outside the region of interest, for examplefn = @(x, y) f(x, y) .* (y <= (1 - x));I = dblquad (fn, 0, 1, 0, 1);​Looking at the Octave source code for dblquad, it seems like it wouldn't be too hard to modify to allow for ymin and ymax being functions of x, which would be a bit more efficient. _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

## Re: help with double quadrature

 In reply to this post by Clinton Winant-2 > On 11 Sep 2017, at 20:17, Clinton Winant <[hidden email]> wrote: > > I have to evaluate an integral over a triangular area: > int_0^1 int_0^y [f(x,y)] dx,dy > > Initially I though I could use dblquad, but is seems the integration limits have to be numbers, not variables.  I have backed of to using quadgk to do the inner integral, then summing over all y.  Does anyone know a better way? Hi, The simplest approach is to use a (singular) affine trasformation that maps your triangle into a square and a quadrature rule that does not place a quadrature node in the point of singularity, e.g. the method descibed in the section "Tensor product-type Gaussian quadrature - Simple but less efficient" in the document at this link:   http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDFTo implement this in Octave you can simply apply dblquad to the transformed integral that you can see in the last formula on page 5 of the above reference. But I think you should switch the backend of dblquad to something that does not place quadrature nodes on the boundary of the domain, e.g. quadgk, rather than quadcc. HTH, c. _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave