How to treat infinite integrals in Octave

I wonder how to translate this below equation to octave:
integral_(-∞)^∞ integral_(-∞)^∞ ( product_(i=1)^6 θ(a_i x + b_i y + c_i)) dx dy
where “θ()” is Heaviside step function or Unit step function.
Also not sure if it requires some manipulations to make it simpler before the integration step…

There are several issues here.

Integrals with infinite bounds are indeed possible. Both quadcc and quadgk handle them for 1-D integrations. So no problems there.

However, the output of the Heaviside step function (Heaviside step function - Wikipedia) is either 0 or 1. If the integrand is finite (i.e., 1) then integrating over an infinite area will always yield an infinite result (Inf) so no reason to even bother with Octave.

If the multiplication of several Heaviside functions results in only a small region of the XY-plane where the integrand is non-zero, then you no longer have an undefined integral. Instead, you have a definite integral and merely need to find x_low, x_high, y_low, y_high which describe the region where the integrand is non-zero. Plotting the function is one easy way to check where it is non-zero. And you don’t have to be precise about determining these limits because the function value outside this area is zero and won’t contribute to the integral. For example, let’s say that function was non-zero in the range [-pi, pi] for both X and Y. Then using limits of [-4, +4] would be fine. For 2-D definite integrals, see the functions integral2 and quad2d.

So I guess the problem could be to simplify it somehow before even trying to enter the equation to octave.
It is enough to set the integration limits to the field size (where the triangles are located). In case field size is not known setting infinite could work better.
Not sure if that numerical integration method would give accurate results (like with 10 decimal points or more).

I made sample case as:

In WFA it seemed to work with below command:
Integrate (Integrate (θ(-1.00x+0.04y+0.29)θ(+0.54x+0.84y+0.29)θ(+0.46x-0.89y+0.29)θ(-1.30x-0.75y+1.66)θ(+0.00x+1.50y+0.09)θ(+1.30x-0.75*y+0.20)) dx x=-0.68…+1.41) dy y=-0.61…+1.34
And the result should be the “overlap” area of the 2 triangles as in below image:

Shifting the other triangle further away should give zero overlap area:
Integrate (Integrate (θ(-1.00x+0.04y+0.29)θ(+0.54x+0.84y+0.29)θ(+0.46x-0.89y+0.29)θ(-1.30x-0.75y+3.20)θ(+0.00x+1.50y-1.04)θ(+1.30x-0.75*y-0.21)) dx x=-0.68…+2.16) dy y=-0.61…+2.09


Attached is an example of running the first integral in Octave. I used anonymous functions to set things up, then plotted the integrand function to get an idea of the limits of integration. The integral2 function throws a warning, because the integrand has discontinuities, but the result is accurate.

H = @(x) ifelse (x > 0, 1, 0) ;
f1 = @(x,y) H(-1.00*x+0.04*y+0.29) ;
f2 = @(x,y) H(+0.54*x+0.84*y+0.29) ;
f3 = @(x,y) H(+0.46*x-0.89*y+0.29) ;
f4 = @(x,y) H(-1.30*x-0.75*y+1.66) ; 
f5 = @(x,y) H(+0.00*x+1.50*y+0.09) ;
f6 = @(x,y) H(+1.30*x-0.75*y+0.20) ;

fint = @(x,y) f1(x,y).*f2(x,y).*f3(x,y).*f4(x,y).*f5(x,y).*f6(x,y) ;

[XX,YY] = meshgrid (-4:0.1:4);
ZZ = fint (XX,YY);
pcolor (XX,YY,ZZ);

[Q, ERR] = integral2 (fint, -1, +1, -1, +1)

tst_hvyside.m (442 Bytes)

The pcolor plot used to find approximate integration limits:

The output from integral2:

warning: quad2d: Maximum number of sub-tiles reached without convergence
warning: called from
    quad2d at line 336 column 7
    integral2 at line 237 column 14
    tst_hvyside at line 15 column 10

Q = 0.1675
ERR = 4.5854e-05


ac =

-0.99912 0.04188 0.28868
0.53583 0.84433 0.28868
0.46330 -0.88620 0.28868
-1.29904 -0.75000 0.54657
0.00000 1.50000 0.36827
1.29904 -0.75000 1.03371
fint=@(x,y) u(ac(1,1)*x+ac(1,2)*y+ac(1,3)).*u(ac(2,1)*x+ac(2,2)*y+ac(2,3)).*u(ac(3,1)*x+ac(3,2)*y+ac(3,3)).*u(ac(4,1)*x+ac(4,2)*y+ac(4,3)).*u(ac(5,1)*x+ac(5,2)*y+ac(5,3)).*u(ac(6,1)*x+ac(6,2)*y+ac(6,3));

warning: quad2d: Maximum number of sub-tiles reached without convergence
warning: called from
quad2d at line 329 column 7
integral2 at line 230 column 14
rot3 at line 373 column 10
q = 0.35814
err = 0.000033944

Looks quite slow though, but answer is given.
octave:2> tic; [q,err]=integral2(fint,bb(1),bb(2),bb(3),bb(4)); toc
warning: quad2d: Maximum number of sub-tiles reached without convergence
warning: called from
quad2d at line 329 column 7
integral2 at line 230 column 14
Elapsed time is 12.4726 seconds.

If the area equation can be converted somehow from integral form to polynomial form that could make it faster.

Always better if you can simplify the problem before going to the computer for calculation. But, these days programmer time is much more valuable than computer time. If a human would take 30 minutes to optimize the calculation while the computer would take 1 minute to make the calculation then the choice is simple.

This would be most likely correct answer for most of the case that computer time is not problem.
But if the calculation is repeating like million times in a loop, that costs could become same, etc.

Of course, but that’s rarely the way scientists interact with Octave. They need to calculate a particular integral in order to move forward with their work. It’s a single calculation, or at worst a few 10’s of runs as various values are explored.

If performance is an issue then an interpreted language like Octave is not appropriate. Better to move to a compiled language and do everything possible to simplify the problem.

If there was some means to simplify the equation, but that seems too hard to be made by myself.
I’m not sure about the scientific calculations, but some case of course it may need only few iterations (rounds) to get it done.
I guess it would then need license for symbolic calculation that could simplify the equations… Not even sure if they can make it or not (like WolframAlpha).
There is really not happening anything to this equation if I try simplify it:
simplify product_(i=1)^6 θ(a_i x + b_i y + c_i)