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

7 次查看（过去 30 天）

Ryan 2023-10-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2039486-integral3-warning-reached-the-maximum-number-of-function-evaluations-10000-the-result-fails-the

编辑： Walter Roberson 2023-11-3

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

9 个评论
显示 7更早的评论隐藏 7更早的评论

Walter Roberson 2023-11-2

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

Ryan 2023-11-2

编辑：Ryan 2023-11-2

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)

请先登录，再进行评论。

请先登录，再回答此问题。