Integral3 Warning: "Reached the maximum number of function evaluations (10000). The result fails the global error test." for a fairly complicated integrand.

Hello,

I've been trying to evaluate a 3D integral at different points on a 2d grid but I appear to run into a lot of issues where the integrand fails and yields a NaN value for seemingly random points. This is strongly affected by various parameters, sometimes the integrals evaluate in some regimes but when I turn some parameters up it fails, even on very lenient relative tolerance values. The integrand is made up of a distribution "cd" (in this case Gaussian), a kernel "kdrift" and a heaviside function "LIMi" that essentially carves out part of the integration domain for tau>0.

u = 10^6; % unit conversion (m to um)
p = 10^15; % unit conversion (s to fs)

%physical constants in SI units 
c = 3*10^8*(u/p); %speed of light (um/ps)

% tau: evaluation time 
tau = 10;
ti = 0;
% energy and velocity
gmi = 10;
bi = (1-1/gmi^2)^(1/2);

% grid parameters
% integral is solved for each point on the grid
% dx is grid scale
% sz defines the number of cells along the axes

dx = 1;
dt = 1;
szt = 10;
sz = 24;

sc = -sz/2*dx:dx:sz/2*dx;
xc = -sz/2*dx:dx:sz/2*dx;

% integral kernels, p denotes primed variable
% kdrift is ordinary longitudinal coulomb field of electron with
% lorentz factor gmi

kdrift = @(x,y,s,xp,yp,sp) (s-sp)./(gmi*(gmi^2*(s-sp).^2+...
    (x-xp).^(2)+(y-yp).^(2)).^(3/2));

% causal limits. Heaviside function that confines fields to the regions in
% which they are causal over the integration
z0i = @(t) c*bi*t;

LIMi = @(x,y,s,xp,yp,sp,t) heaviside(-c^2*(t-ti).^2+(x-xp).^2+(y-yp).^2+(s-sp+z0i(t)).^2);

% charge distribution (Gaussian)
% bunch dimensions coupled to grid scale dx (um)

sigX = 5*dx;
sigY = sigX;
sigZ = sigX;

N = 10^9;
    
cd = @(xp,yp,sp) (N/((2*pi)^(3/2)*sigX*sigY*sigZ))...
    *exp(-xp.^2/(2*(sigX)^(2)))...
    .*exp(-yp.^2/(2*(sigY)^(2)))...
    .*exp(-sp.^2/(2*(sigZ)^(2)));

% 2D grid

Efieldz = zeros(sz+1,sz+1);

% integration limits
DX = 5*sigX;
DY = 5*sigY;
DZ = 5*sigZ;

% loop solves integral for each point in the grid

for i = 1:1:sz+1
    for k = 1:1:sz+1

        Efieldz(i,k) = integral3(@(xp,yp,sp) cd(xp,yp,sp)...
        .*kdrift(xc(i),0,sc(k),xp,yp,sp)...
        .*LIMi(xc(i),0,sc(k),xp,yp,sp,tau)...
        ,-DX,DX,-DY,DY,-DZ,DZ,'RelTol',1e-2,'method','tiled');   
               
    end      
end

f1 = figure;
imagesc(sc,xc,Efieldz,[min(Efieldz,[],'all'), max(Efieldz,[],'all')]),colorbar;
grid on;
title('EFieldz');
xlabel('s');
ylabel('x');
hold on
axis square;

The warning I get is

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

> In integral2Calc>integral2t (line 129)

In integral2Calc (line 9)

In integral3>innerintegral (line 137)

In integral3>@(x)innerintegral(x,fun,yminx,ymaxx,zminxy,zmaxxy,integral2options) (line 111)

In integralCalc/iterateScalarValued (line 323)

In integralCalc/vadapt (line 132)

In integralCalc (line 75)

In integral3 (line 113)

which appears to correspond to the points at which the integrations fail. The integrations work for tau=0 and gmi=10, but fail for higher gmi (such as 100) or tau=1 and above. Tau being non zero essentially "switches on" the heaviside function otherwise it should be 1 across the whole domain for tau=0. I've tried different coordinate systems such as spherical or cynlindrical which appear to help slightly but I still get broadly the same warnings and NaN results. I'm really not sure why I'm having these problems, I can see the kernel going singular at some points but it isn't the nastiest integration. It's confusing to me that it succeeds at some points but fails at others, especially as these points have nothing distinct about them. My colleague tried to run it in mathematica and had even more issues there!

EDIT:

Amended with plotting function, for tau<ti and gmi=1 (gmi cannot <1) I obtain

[Walter Roberson]

When tau=1 you have discontinuity and need to use iterated instead of tiled

[Ryan]

Thanks for the comment Walter. It appears to behave slightly differently on iterated but I'm still having more or less the same problem!

[Walter Roberson]

The code does not define ti

[Walter Roberson]

I defined ti arbitrarily, and then I modified the code to create symbolic numbers and to use vpaintegral() . Even just the first of the entries is taking a fair bit of time.

I tried putting it through Maple, but unfortunately I encounted a bug in Maple when I did that (which I have now reported to Maplesoft)

[Ryan]

Ah I must have removed it, ti should equal zero. Generally tau>=ti

I was cleaning up the code to put it into here because there was a lot of clutter, I must have accidentally put ti on the chopping block! I have edited the code. Do you think VPAintegral might be better suited to this problem numerically?

[Ryan]

Going to give this a cautious bump as I have yet to resolve these issues

[Walter Roberson]

I was able to use vpaintegral for a few locations but it was slow enough that I gave up.

The Maple developers told me yesterday that they have fixed the integration bug I encountered in Maple when trying this, but it will probably be December before that version is available for testing.

[Walter Roberson]

Could you indicate a couple of specific locations that it fails at?

[Ryan]

Hi Walter, thanks for your further comments. I've amended my post with plotting script and an example of how the integral fails in a limiting case. Here tau<ti so my heaviside LIMi should be equal to one over the domain. I've also shifted where I evaluate my integral by adjusting sz to 24. gmi=1 means in the denominator of my kernal I have (s-sp)^2+(x-xp)^2+(y-yp)^2. In this case there is a bit of a pattern to the failure, it's all down the axis s=0 which could offer a clue. I'm wondering if perhaps it's something to do with integral3 struggling in the vicinity of the singularity due to the steepness of the gradient? Not the singularity itself persay but precision problems relating to large numbers sampled close to it? (I'm not terribly familiar with the inner workings of integral3)