Why does my code get stuck when trying to convert a double integral into a double?

syms x y z

f(x,y) = sqrt(1+((4*y)/(121*sqrt((4*x^2)/121 + (4*y^2)/121))^2)+((4*x)/(121*sqrt((4*x^2)/121 + (4*y^2)/121))^2))
minX= -sqrt(220.523)
maxX= sqrt(220.523)   %input("Limite superior de X:")
intX = int(f,x,minX,maxX)
minY= -sqrt(220.523-x^2)  %input("Limite inferior de Y:")
maxY= sqrt(220.523-x^2)   %input("Limite superior de Y:")
intY = int(f,y,minY,maxY)
integralDobleTipoUno = int(intY,x,minX,maxX)
integralDobleTipoDos= int(intX,y,minY,maxY)
var = double(integralDobleTipoUno)

When reaching "var", it loads forever

Answer 1

[Walter Roberson]

if you are after the numeric solution then switch to vpaintegral()

[Dyuman Joshi]

syms x y z

f(x,y) = sqrt(1+((4*y)/(121*sqrt((4*x^2)/121 + (4*y^2)/121))^2)+((4*x)/(121*sqrt((4*x^2)/121 + (4*y^2)/121))^2));
minX= -sqrt(220.523);
maxX= sqrt(220.523);   %input("Limite superior de X:")
intX = int(f,x,minX,maxX);
minY= -sqrt(220.523-x^2);  %input("Limite inferior de Y:")
maxY= sqrt(220.523-x^2);   %input("Limite superior de Y:")
intY = int(f,y,minY,maxY);
integralDobleTipoUno = vpaintegral(intY,x,minX,maxX)

integralDobleTipoUno = 

integralDobleTipoDos= int(intX,y,minY,maxY)

integralDobleTipoDos = 

%var = double(integralDobleTipoUno)

However, you won't be able to do the same for the 2nd integral, because it's output is in terms of symbolic variables, not symbolic numerical values.

[Walter Roberson]

@Dyuman Joshi

You can keep the call to double(), because it is operating on integralDobleTipoUno which vpaintegral() has already reduced to numeric form

format long g
syms x y z

f(x,y) = sqrt(1+((4*y)/(121*sqrt((4*x^2)/121 + (4*y^2)/121))^2)+((4*x)/(121*sqrt((4*x^2)/121 + (4*y^2)/121))^2));
minX= -sqrt(220.523);
maxX= sqrt(220.523);   %input("Limite superior de X:")
intX = int(f,x,minX,maxX);
minY= -sqrt(220.523-x^2);  %input("Limite inferior de Y:")
maxY= sqrt(220.523-x^2);   %input("Limite superior de Y:")
intY = int(f,y,minY,maxY);
integralDobleTipoUno = vpaintegral(intY,x,minX,maxX)

integralDobleTipoUno = 

integralDobleTipoDos= vpaintegral(intX,y,minY,maxY)

integralDobleTipoDos = 

var = double(integralDobleTipoUno)

var = 
           692.792994216759 +  8.32119535172846e-05i

[Walter Roberson]

integralDobleTipoUno and IntegralDobleTipoDos are intended to be the same overall integrals, just with the integration order reversed -- one integrates by x first and the ny, and the other integrates by y first and then x.

But the boundary conditions for y depend on x, so if you integrate by x first there is no remaining unbound x to be able to carry out the integration on y. That is why integralDobleTipoDos is not able to find a result. Integration order matters when the boundaries are dependent on a variable of integration.

[Guillermo]

@Walter Roberson

What I don't understand is how apps like Wolfram alpha give me numeric answers with the same integrals (i.e. for this same double integral the result in Wolfram is 704.151), while MATLAB either gets stuck calculating the double, or the vpaintegral has an i in it.

[Walter Roberson]

For one thing, remember GIGO, "garbage in, garbage out".

When you write 220.523 you are going to get some number between approximately 220523/1000 ± 10^-12. Not necessarily exactly the same number for Windows as for Mac Intel or Mac ARM or Linux Intel. I think it should be consistent between Windows Intel and Windows AMD but I would not want to guarantee.

Anyhow you get some particular binary double precision floating point number. You do not get 220523/1000.

And then you sqrt() that as a double precision floating point number. The result can be different for Intel vs AMD vs ARM, and in theory could be different for different models within those architectures, or could in theory depend on whether it was stand-alone or part of a vectorized operation.

Then you apply those increasingly uncertain numbers as the boundaries of integration of an integral that has a division, and you wonder why different software systems give you different answers.

int() is a request to produce indefinitely precise results, exact values, algebraic numbers if possible, but you are using inexact boundaries. You were never going to get back meaningful results.

[Walter Roberson]

format long g
syms x y z

reltol = 1e-10;

magic_number = sym(220523)/1000;
f(x,y) = sqrt(1+((4*y)/(121*sqrt((4*x^2)/121 + (4*y^2)/121))^2)+((4*x)/(121*sqrt((4*x^2)/121 + (4*y^2)/121))^2));
maxX = sqrt(magic_number);   %input("Limite superior de X:")
minX = -maxX;
intX = int(f,x,minX,maxX);
maxY = sqrt(magic_number-x^2);   %input("Limite superior de Y:")
minY = -maxY;
intY = int(f,y,minY,maxY);
integralDobleTipoUno = vpaintegral(intY, x, minX, maxX, 'RelTol', reltol)
integralDobleTipoDos= vpaintegral(intX, y, minY, maxY, 'reltol', reltol)
var = double(integralDobleTipoUno)

gives

692.79299543899 + 8.42632681212178e-05i

If you take form of int(intY, x, minX, maxX) and ask Wolfram Alpha for the value

int(int(((4*x)/(484*x^2 + 484*y^2) + (4*y)/(484*x^2 + 484*y^2) + 1)^(1/2), y, -(220523/1000 - x^2)^(1/2), (220523/1000 - x^2)^(1/2)), x, -(10^(1/2)*220523^(1/2))/100, (10^(1/2)*220523^(1/2))/100)

then it gives 692.793 + 0.0000938591

[Guillermo]

@Walter Roberson

Crazy to think evven my CPU can make a difference, although it makes sense. Anyhow, thank you so much!

[Walter Roberson]

format long g
syms x y z

magic_number = sym(220523)/1000;
f(x,y) = sqrt(1+((4*y)/(121*sqrt((4*x^2)/121 + (4*y^2)/121))^2)+((4*x)/(121*sqrt((4*x^2)/121 + (4*y^2)/121))^2));
maxX = sqrt(magic_number);   %input("Limite superior de X:")
minX = -maxX;
intX = int(f,x,minX,maxX);
maxY = sqrt(magic_number-x^2);   %input("Limite superior de Y:")
minY = -maxY;
intY = int(f,y,minY,maxY);
Sf = simplify(f, 'steps', 50)

Sf(x, y) = 

Tx = 0; Ty = -1/200;
combine(subs(maxY, x, Tx))

ans = 

That tells us that for x = 0, that y = -1/200 is well within range

temp = simplify(limit(f, x, Tx), 'steps', 50)

temp(y) = 

y_bound = solve(children(temp, 1))

y_bound = 

limit(limit(f, x, Tx), y, Ty)

ans = 

but we got a complex-valued result.

This tells us that the complex component we got was not an accident: there is a range of values where f goes complex.