Double symbolic integral - "Explicit integral could not be found"

At i=1 and j=3 the following is displayed "Warning: Explicit integral could not be found." A preliminary evaluation revealed that F(i=1,j=4)is failing to integrate x over the specified range with the result still containing the variable x.

Results:

F(i=1,j=2) "No Problem"
F = -(154946195819313285227594878996707*y*exp(-y^2)*((7186705221432913*2^(1/2)*pi^(1/2)*(erf((2^(1/2)*y)/2) + 1))/36028797018963968 - 1)^2)/81129638414606681695789005144064

F(i=1,j=3) "No Problem"
F = (1113552634535825382050544662406231220558348417491*pi^(1/2)*y*exp(-y^2/2)*(erf(y) + 1)*(7186705221432913*2^(1/2)*pi^(1/2) + 7186705221432913*2^(1/2)*pi^(1/2)*erf((2^(1/2)*y)/2) - 36028797018963968))/52656145834278593348959013841835216159447547700274555627155488768

F(i=1,j=4) "Warning: Explicit integral could not be found." 
F = int((154946195819313285227594878996707*x*y*exp(-x^2/2)*exp(-y^2/2)*((7186705221432913*2^(1/2)*pi^(1/2)*(erf((2^(1/2)*x)/2) + 1))/36028797018963968 - (7186705221432913*2^(1/2)*pi^(1/2)*(erf((2^(1/2)*y)/2) + 1))/36028797018963968)^2)/81129638414606681695789005144064, x == -Inf..y)

Code:

clear all
clc

n = 4;
syms i j x y t w;

C2 = (2*pi())^(-0.5);

C3 = factorial(n)/(factorial(i-1)*factorial(j-i-1)*factorial(n-j));

fx = C2*exp(-0.5*x^2);
fy = C2*exp(-0.5*y^2);

ft = C2*exp(-0.5*t^2);
fw = C2*exp(-0.5*w^2);

Fx = int(ft,-inf,x);
Fy = int(fw,-inf,y);

for i=1:n/2
    for j=(i+1):3
        EC3 = eval(C3);
        Fxy = EC3*x*y*fx*fy*Fx^(i-1)*(Fy-Fx)^(j-i-1)*(1-Fy)^(n-j);
        F   = int(Fxy,x,-inf,y);
        Exy(i,j) = double(int(F,y,-inf,inf)) 
    end
end

[Walter Roberson]

Do you have reason to believe that there is a closed form integral for that?

[Nick Martin]

Given that all of the components are continuous on the interval, yes, I do believe that there is a closed form integral. It's a standard equation for calculating the covariance of two Gaussian random variables.

[Walter Roberson]

Continuity is a really insufficient criteria for closed form integrals to exist. For example an ellipse is continuous but extensive study over centuries has found no closed form integral for its arc length. http://en.m.wikipedia.org/wiki/Arc_length

Maple for one finds no closed form integral for the case you raise here.

[Nick Martin]

Good point and thanks for evaluating with Maple. I'm going to investigate an alternative method to evaluating/approximating this integral.

Answer 1

[Walter Roberson]

If that integral has a closed form at all, then it is difficult to find. MATLAB warns you that it has left an integral unevaluated. But it is only a warning. When it comes time to do the double() on the nested unevaluated integrals, a numeric integration will be done to approximate the value. The warning helps you to be aware that the numeric result you get out will be an approximation whereas the other results are exact to within floating point round-off.

If you really do not want to see the warning, then you can use woff(). See lastwarning() to find the identifier of the warning to be able to turn it off.

[Nick Martin]

I'm comfortable with an approximated result, but the problem is that the code is interrupted due to the presence of x. The next step is evaluating the integral over (dy) and the presence of x locks it up.

[Walter Roberson]

The unevaluated symbolic integral, integrated, will result in a nested unevaluated symbolic integral. And that's OK because double() tells matlab to find the numeric result.

It is not an error for an integral to be symbolically unevaluated. For example

int(int(f(x, y), x), y)

Is a perfectly valid symbolic construct and can be manipulated symbolically. The symbolic engine is happy to keep it in the form of an unevaluated integral.

[Nick Martin]

Walter - I appreciate the feedback, but unfortunately MATLAB is failing to evaluate the integral. I tried using double (multiple errors) and vpa (result was numeric::int ?) - see below. Any advice on what I need to change to obtain a valid result?

F = int(Fxy,x,-inf,y) Warning: Explicit integral could not be found.

F =

int((6*x*y*exp(-0.5*x^2)*exp(-0.5*y^2)*(0.5*erf(0.70710678118654752440084436210485*x) - 0.5*erf(0.70710678118654752440084436210485*y))^2)/pi, x == -Inf..y)

Exy = double(int(F,y,-inf,inf));

Warning: Explicit integral could not be found. Warning: Explicit integral could not be found. Error using mupadmex Error in MuPAD command: DOUBLE cannot convert the input expression into a double array. If the input expression contains a symbolic variable, use the VPA function instead.

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

G = int(F,y,-inf,inf);

Exy = vpa(G,5); Warning: Explicit integral could not be found. Exy

Exy =

numeric::int(numeric::int((6*x*y*exp(-0.5*x^2)*exp(-0.5*y^2)*(0.5*erf(0.70711*x) - 0.5*erf(0.70711*y))^2)/pi, x == -Inf..y), y == -Inf..Inf)

[Walter Roberson]

This is potentially a bug in how vpa() does its work. You could try,

G = evalin(symengine, 'numeric::int', F, sym('y=-infinity..infinity'));
Exy = double(G);

but it may use up all of your memory.

[Walter Roberson]

When I ask Maple to do a numeric integration using 13 digits or fewer, it finds a solution in not too long, of about -0.9549296585514. But when I ask it to do 14 or 15 digits, it returns the integral unevaluated. When I ask for 16 digits it runs through all of my memory and demands more.

When I plot I do not see anything about the values that would lead to that behaviour. The peak is at roughly y = -1.2, x = -1.9 and is positive but the values slip negative fairly close-by so it is plausible that the overall integral is negative.