Double symbolic integral - "Explicit integral could not be found"
3 次查看(过去 30 天)
显示 更早的评论
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
4 个评论
Walter Roberson
2015-6-12
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.
回答(1 个)
Walter Roberson
2015-6-12
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.
5 个评论
Walter Roberson
2015-6-14
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
2015-6-14
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.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!