Search code examples
matlabcalculusnumerical-integration

Errors when using the Integral2 function in MATLAB


As far as I can tell, no one has asked this. I've been asked to compute the double integral of a function, and also the same double integral but with the order of integration swapped (i.e: first integrate for dydx, then dxdy). Here is my code:

    %Define function to be integrated
f = @(x,y) y^2*cos(x);

    %First case. Integration order: dydx
ymin = @(x) cos(x);
I = integral2(f,ymin,1,0,2*pi)

    %Second case. Integration order: dxdy
xmin = @(y) asin(y)+2*pi/2;
xmax = @(y) asin(y)-pi/2;
B = integral2(f,xmin,xmax,-1,1)

The error I'm getting is this:

Error using integral2 (line 71)

XMIN must be a floating point scalar.

Error in EngMathsA1Q1c (line 5)

I = integral2(f,ymin,1,0,2*pi)

I'm sure my mistake is something simple, but I've never used Integral2 before and I'm lost for answers. Thank you.


Solution

  • Per the integral2 documentation, the variable limits are given as the second pair of limits. So your first integral should be

    %    Define function to be integrated
    f = @(x,y) y.^2.*cos(x);
    
    %    First case. Integration order: dydx
    ymin = @(x) cos(x);
    I = integral2(f,0,2*pi,ymin,1);
    

    The set of constant limits always goes first, and Matlab assumes the first argument of f is associated with the first set of limits while the second argument of f is associated with the second set of limits, which may be a function of the first argument.


    I point out that second part because if you wish to switch the order of integration, you also need to switch the order of the inputs of f accordingly. Consider the following example:

    fun = @(x,y) 1./( sqrt(2*x + y) .* (1 + 2*x + y).^2 )
    

    A nice little function that is not symmetric in its arguments (i.e., fun(x,y) ~= fun(y,x)). Let's integrate this over an elongated triangle in the first quadrant with vertices at (0,0), (2,0), and (0,1). Then integrating with dA == dy dx, we have

    >> format('long');
    >> ymax = @(x) 1 - x/2;
    >> q = integral2(fun,0,2,0,ymax)
    q =
       0.220241017339352
    

    Cool. Now let's integrate with dA == dx dy:

    >> xmax = @(y) 2*(1-y);
    >> q = integral2(fun,0,1,0,xmax)
    q =
       0.241956050772765
    

    Oops, that's not equal to the first calculation! That's because fun is defined with x as the first argument and y as the second, but the previous call to integral2 is implying that y is the first argument to fun, and it has constant limits of 0 and 1. How do we fix this? Simply define a new function that flips the arguments:

    >> fun2 = @(y,x) fun(x,y);
    >> q = integral2(fun2,0,1,0,xmax)
    q =
       0.220241017706984
    

    And all's right with the world. (Although you may notice small differences between the two correct answers due to the error tolerances of integral2, which can be adjusted via options per the documentation.)