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 |
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 |
On Mon, Sep 11, 2017 at 4:04 PM, Juan Pablo Carbajal <[hidden email]> wrote:
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 |
On Mon, Sep 11, 2017 at 3:31 PM, Nicholas Jankowski <[hidden email]> wrote:
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 |
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));_______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
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 |
Free forum by Nabble | Edit this page |