Computational Science People:

The title is the question. I am trying to numerically compute a certain integral over a square in Octave (an open-source Matlab clone). I'm getting the wrong answer. To show that something really is going wrong, I've created the following toy problem: Let $f = \chi_{(-\infty, 1)}$ ($f(t) = 1$ if $t < 1$, otherwise $0$). Then $\int_0^3 \int_0^3 f(t) f(t')\,dt'\,dt = 1$. I wrote an Octave script to compute this:

    junk=10;  %this is here just so Octave doesn't complain about the file name
    function out =R1(t)   
    if t < 1     
       out = 1; 
       out =0;    

    dIF = @(t1,t2) ones(size(t1))* R1(t1(1))*R1(t2);

    dblquad(dIF, 0, 3, 0, 3)  %should be 1.  Octave incorrectly gives 3.

    [xx, yy] = meshgrid (tx, ty);

    for i=1:31
      for j=1:31 
      end%for j
    end%for i   
    mesh(xx,yy,z) %graph is correct  

 %end script

Octave returns $3$ for the value of the double integral, and it should be $1$. As you can see, the script plots the integrand, which is, as intended, $1$ if $x$ and $y$ are both between $0$ and $1$, and $0$ otherwise. I defined the function dIF the way I did because in my actual problem, I need to pass additional parameters to the function I'm integrating, and this is the only way I know how to do it.

In my actual problem, I'll probably just cut my double integral into nine pieces, not requiring piecewise functions, and add them. But I'd like to know if there's a better way to do it, in case I face a similar problem in the future.

I have not attempted this in Matlab, only Octave. I tested it just now with a copy of Octave I downloaded today from what I am quite sure is a reputable source.


1 Answer 1


Your R1 function should be vector valued (it appears as though it is scalar valued). For example, this appears to work:

function out = R1(t)   
  out = t < 1;

dIF = @(x,y) R1(x) .* R1(y);

dblquad(dIF, 0, 3, 0, 3)
  • $\begingroup$ I got it to work for my "real" problem, and it kind of worked: the relative error was about 1%. I ended up cutting the domain (a square) into 9 rectangles, doing double integrals of non-piecewise functions, and adding. Due to symmetries of the problem, I only needed three integrand functions and I only needed to compute four double integrals. $\endgroup$ Commented Jul 22, 2013 at 23:45
  • $\begingroup$ Glad to hear that it worked out... $\endgroup$ Commented Jul 24, 2013 at 3:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.