help with double quadrature

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

help with double quadrature

Clinton Winant-2
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?

_______________________________________________
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

Juan Pablo Carbajal-2
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

NJank
On Mon, Sep 11, 2017 at 4:04 PM, Juan Pablo Carbajal <[hidden email]> 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

Dmitri A. Sergatskov
On Mon, Sep 11, 2017 at 3:31 PM, Nicholas Jankowski <[hidden email]> wrote:
On Mon, Sep 11, 2017 at 4:04 PM, Juan Pablo Carbajal <[hidden email]> 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?



​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

Nir Krakauer-3
Matlab's documentation suggests using dblquad with the function modified to be 0 outside the region of interest, for example

fn = @(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

Carlo de Falco-2
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.PDF

To 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