Using symbolic integration: Why does my double integral returns different numerical values depending on integration order?

I have the double integral . Using int and double to numerically evaluate this I get:

   syms x y
   fun1 = @(x, y) (112).*exp(-7.*x-9.*y);
   a = int(int(fun1, y, 1/4), y, 0, 1);
   double(a)
   ans = 0.6911

If I change the order of integration and integrate , I get the correct answer.

   syms x y 
   fun = @(x, y) (112).*exp(-7.*x-9.*y);
   b = int(int(fun, y, 0, x), x, 0, 1/4);
   double(b)
   ans = 0.7053

Why is this the case?

Answer 1

[Paul]

Comment: I'm not going to use the anonymous functions because, as far as I know, those are invalid inputs to int and I don't know how int is working with those inputs.

In short: the limits of integration in the two integrals do not cover the same area in the x-y plane.

Because the question states that that second integral is correct, let's start with that

syms x y
f(x,y) = (112).*exp(-7.*x-9.*y);
I1 = int(int(f(x,y),y,0,x,'hold',true),x,0,sym(1)/sym(4),'hold',true)

I1 = 

I1 = release(I1)

I1 = 

double(I1)

ans = 0.7053

I1 is the integral of f(x,y) over the blue triangular area

area(0:.01:.25,0:.01:.25)

So, if we want reverse the order of integration then we see that for each value of y, x varies from y to 1/4, and y varies from 0 to 1/4. However, the question shows the limits of integration on y from 0 to 1.

I2 = int(int(f(x,y),x,y,sym(1)/sym(4),'hold',true),y,0,sym(1)/sym(4),'hold',true)

I2 = 

I2 = release(I2)

I2 = 

simplify(I1 - I2)

ans = 

0

Answer 2

[Walter Roberson]

   syms x y
   fun1 = @(x, y) (112).*exp(-7.*x-9.*y);
  int(fun1, y, 1/4)

ans = 

Which variable did it integrate with respect to? Answer: x by default, since x is the first variable returned by symvar()

  symvar(sym(fun1))

ans = 

So we can put the variable in explicitly

a = int(int(fun1, x, y, 1/4, 'hold', true), y, 0, 1, 'hold', true)

a = 

and that matches the formula you asked to integrate.

Now let us try your other version:

   fun = @(x, y) (112).*exp(-7.*x-9.*y);
   b = int(int(fun, y, 0, x, 'hold', true), x, 0, 1/4, 'hold',true)

b = 

Is that the same? Notice that the -7 here matches to the integration range 0 to 1/4, whereas the original that you said you wanted to integrate, the -7 is associated with the integration range y to 1/4 .

The two are not equivalent.

[Jame Tran]

Hi Walter,

I still get the same wrong numeric answer if I do.

fun1 = @(x, y) (112).*exp(-7.*x-9.*y);
a = int(int(fun1, x, y, 1/4), y, 0, 1);
double(a)
ans =  0.6911

The inner intergral is with respect to x, so originally integrating using default variable x makes sense. If I take the same integral and swap the integration order to dydx, I get 0.7053 which is the correct answer for the integral.

[Walter Roberson]

syms x y
fun1 = @(x, y) (112).*exp(-7.*x-9.*y)

fun1 = function_handle with value:
    @(x,y)(112).*exp(-7.*x-9.*y)

a = int(int(fun1, x, y, 1/4, 'hold', true), y, 0, 1, 'hold', true)

a = 

b = int(int(fun1, y, 0, 1, 'hold', true), x, y, 1/4, 'hold', true)

b = 

ar = simplify(release(a))

ar = 

br = simplify(release(b))

br = 

double(ar)

ans = 0.6911

double(br)

Error using symengine
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for variables.

Error in sym/double (line 800)
Xstr = mupadmex('symobj::double', S.s, 0);

What happened? Well, when you swap the order of integration without rewriting the expression, then y in the limit for the integral x becomes a "different" y then the y being integrated over first.

[Jame Tran]

Hi Walter,

I'm not quite understanding what you're saying. I'm changing the order of integration, but I also rewrote the the bounds when changing the order. Both integrals are over the triangular area defined by y < x < 1/4, 0 < y < 1 right?

[Walter Roberson]

In the first version, fun1, the dependent variable whose bounds are expressed in terms of the other variable, is associated with the -7 multiplier in the exp()

In the second version, fun, the dependent variable whose bounds are expressed in terms of the other variable, is associated with the -9 multiplier in the exp()

Unless you are saying that you believe that they should calculate the same result through two different methods? I would have to think about that a bit.

[Jame Tran]

Hi Walter,

I am saying they should calculate the same result, since they're integrating the same function over the same triangular area right? I'm still a student so I might be misunderstanding something haha.

[Walter Roberson]

Both integrals are over the triangular area defined by y < x < 1/4, 0 < y < 1 right?

No. Your second integral implies y <= x, but your first integral implies x <= y