1
$\begingroup$

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; 
    else  
       out =0;    
    end%if                               
    end%function

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

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

    tx=0:0.1:3;
    ty=0:0.1:3;
    [xx, yy] = meshgrid (tx, ty);
    z=zeros(31,31);

    for i=1:31
      for j=1:31 
        z(i,j)=dIF(xx(i,j),yy(i,j)); 
      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.

$\endgroup$
2
$\begingroup$

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;
end

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

dblquad(dIF, 0, 3, 0, 3)
$\endgroup$
  • $\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$ – Stefan Smith Jul 22 '13 at 23:45
  • $\begingroup$ Glad to hear that it worked out... $\endgroup$ – Matthew Emmett Jul 24 '13 at 3:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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