Numerical integration and processing time

Hello,

I need to evaluate the following integral A. So, I devided it into three functions as

The functions of interest are built as

function [out] = funct1(x,y)
ny=length(y);
out=zeros(1,ny);
for i=1:ny
    S1=@(v,theta) v.*exp(-v).*exp(-sqrt(v.^2+x.^2-2.*x.*v.*cos(theta)));
    S2=@(theta) exp(-sqrt(y(i).^2+x.^2-2.*x.*y(i).*cos(theta)));
    out(i)=integral(S2,0,2*pi,'ArrayValued',true).*exp(-quad2d(S1,0,y(i),0,2*pi));
end
end

function [out] = funct2(x)
nx=length(x);
out=zeros(1,nx);
for i=1:nx
    S3=@(y) y.*exp(-y).*funct1(x(i),y);
    out(i)=integral(S3,0,inf,'ArrayValued',true);
end
end

function [out] = funct3(a)
S=@(x) funct2(x);
out=a*integral(S,0,inf);
end

One major problem is that such codes give a very delayed result with warning message as

tic
A=funct3(1)
toc
%%%%%%%%%%%%%%%%%%%%%
Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.
A =

    2.3982

Elapsed time is 75.049599 seconds.

Any help please to reduce such processing delay and avoid the warning message.

Many thanks in advance.

[Walter Roberson]

The formula for the second integral is not correct.

The second-last integral is with respect to dx, and does not use x as a limit. Therefore after the integration, there will be no remaining x in the equation. But the outer integral uses x, and the outer integral is with respect to dy, so you would have unbound x in the result.

It would make more sense of the order of integration was dy dx rather than dx dy.

[Yassine Hmamouche]

thanks, it's corrected!

[Walter Roberson]

I have been working on this using Maple. I find that if the number of digits of precision is set too low then Maple gives up after a while and returns the expression unevaluated, which is the way it signals that it detected that the integration error was too high, that the integral did not convert numerically for those parameters at that precision. As I increase the precision, it starts taking a long long time.

For exploration I found it helped to break up the code:

syms x y v theta
assume(x>=0 & y>=0 & v >= 0 & theta >=0)
p45 = int(int(v*exp(-v)*exp(-sqrt(x^2 + v^2 - 2*x*v*cos(theta))), theta, 0, 2*pi), v, 0, y);
p3 = int(exp(-sqrt(x^2 + y^2 - 2*x*y*cos(theta))), theta, 0, 2*pi);

p2 = int(y*exp(-y)*p3*exp(-p45), y, 0, inf);
p1 = int(x*exp(-x)*p2, x, 0, inf);

You will want to execute a line at a time, as the p2 calculation runs indefinitely. When I explore in Maple, I can use Int() instead of int() : Int() is an unevaluated integration that can be further manipulated symbolically.

Unfortunately so far I have not been able to find any changes that would obviously made a difference.

I am running into theoretical issues in analyzing at infinity. For example,

sqrt(x^2 + y^2 - 2*x*y*cos(theta))

As x -> infinity, for positive y and cos(theta) > 0, that gives an expression that is approximately infinity^2 - infinity which is undefined. You have to do a limit process (which would say that infinity^2 increases much faster than 2*infinity*y*cos(theta) decreases, so as a limit you have like sqrt(infinity^2) . Unfortunately the limit() functions do not like to process the entire expression...

[Yassine Hmamouche]

I considered relative error in integrals as "RelTol",1e-4 and reduced the interval of integration from inf to 1000, I get

tic
A=funct3(1)
toc

A =

    2.3982

Elapsed time is 38.308114 seconds.

Which eliminates the warning message and reduces the previous processing time by half.

I wonder if there is any way to vectorize funct 1, 2 and 3, which will further shrink the processing time.

Thanks in advance!

[Walter Roberson]

For example for x = 10^5, then p4 works out as 6.6951418857737784647*10^(-43424)... after several hours.

[Walter Roberson]

For the overall integral, coding in terms of vpaintegral() with RelTol 1e-4, then MATLAB comes up with 0.5485 after 14307 seconds (just under 4 hours)

The code was

syms x y v theta
assume(x>=0 & y>=0 & v >= 0 & theta >=0)
Int = @(F,var,lb,ub) vpaintegral(F, var, lb, ub, 'RelTol',1e-4);
p5 = Int(v*exp(-v)*exp(-sqrt(x^2 + v^2 - 2*x*v*cos(theta))), theta, 0, 2*pi);
p45 = Int(p5, v, 0, y);
p3 = Int(exp(-sqrt(x^2 + y^2 - 2*x*y*cos(theta))), theta, 0, 2*pi);

tic; p2 = Int(y*exp(-y)*p3*exp(-p45), y, 0, inf); toc
tic; p1 = Int(x*exp(-x)*p2, x, 0, inf); toc     %4 hours !!!
digits(16);
tic; funct3 = vpa(p1); toc
disp(funct3)