A double symbolic integral in Matlab

Question

I am trying to calculate some integrals, for example:

a = 1/sqrt(2); 
b = -5;
c = 62;
d = 1;
f = exp(-x^2-y^2)*(erfc((sym(a) + 1/(x^2+y^2)*(sym(b)*x+sym(d)*y))*sqrt((x^2+y^2)*sym(10.^(c/10))))...
                 + erfc((sym(a) - 1/(x^2+y^2)*(sym(b)*x+sym(d)*y))*sqrt((x^2+y^2)*sym(10.^(c/10)))));    
h = int(int(f,x,-Inf,Inf),y,-Inf,Inf);

It will occur error like this:

Warning: Explicit integral could not be found.

Then, I try to use vpato calculate that integral,and get the result like this

vpa(int(int(f,x,-Inf,Inf),y,-Inf,Inf),5)
numeric::int(numeric::int(exp(- x^2 - y^2)*(erfc(((6807064429273519*x^2)/4294967296 + (6807064429273519*y^2)/4294967296)^(1/2)*(2^(1/2)/2 + (5*x - y)/(x^2 + y^2))) + erfc(((6807064429273519*x^2)/4294967296 + (6807064429273519*y^2)/4294967296)^(1/2)*(2^(1/2)/2 - (5*x - y)/(x^2 + y^2)))), x == -Inf..Inf), y == -Inf..Inf)

I already tried to change the interval [-Inf,Inf] to [-100,100], and get the same above result:

numeric::int(numeric::int(exp(- x^2 - y^2)*(erfc(((6807064429273519*x^2)/4294967296 + (6807064429273519*y^2)/4294967296)^(1/2)*(2^(1/2)/2 + (5*x - y)/(x^2 + y^2))) + erfc(((6807064429273519*x^2)/4294967296 + (6807064429273519*y^2)/4294967296)^(1/2)*(2^(1/2)/2 - (5*x - y)/(x^2 + y^2)))), x == -100..100), y == -100..100)

My question is why vpa in this case could not return to a real value? There are something wrong in above Matlab code? (I, myself, could not find the bug so far) Thank you in advance for your help.

horchler horchler · Accepted Answer · 2016-04-11T17:49:14

It is unlikely that there is an analytic solution to this integral so using int may not be a good choice. In some cases vpa can be used to for a numeric solution. When this fails (by returning a call to itself) it may be for several reason: the integral may not exist, the integral may converge too slowly, singularities may cause issues, the integrand may be highly oscillatory or non-smooth, etc. Mathematica also struggles with this integral.

You can try calculating the integral numerically using integral2:

a = 1/sqrt(2);
b = -5;
c = 62;
d = 1;
f = @(x,y)exp(-x.^2-y.^2).*(erfc((a + 1./(x.^2+y.^2).*(b*x+d*y)).*sqrt((x.^2+y.^2)*10^(c/10)))...
          +erfc((a - 1./(x.^2+y.^2).*(b*x+d*y)).*sqrt((x.^2+y.^2)*10^(c/10))));
h = integral2(f,-Inf,Inf,-Inf,Inf)

which returns 5.790631184403967. This compares well with Mathematica's numerical integration using NIntegrate. You can try specifying smaller absolute and relative tolerances for integral2 to get more accurate values, but this will result in much slower compute times.

A double symbolic integral in Matlab

1 Answers