Why does this integral of a relatively well behaved function return NaN values?

Hello,

I am trying to evaluate an integral of a function (as part of a cross-check between my own results and the results from a piece of literature), however I'm struggling to get MATLAB to integrate what for all intents and purposes is a relatively well behaved function (though admittedly it is a fairly complicated one).

See the example code below which I have streamlined a bit for readability. I setup some initial parameters (physical dimensions and a dimensionless energy or lorentz factor), setup a positional array and then rescale it as well as defining a dimensionless parameter, eta. Then I simply run the integral for different values of the dimensionless position array d. (I'm aware I could use the array valued feature of integral to do it without the loop but it's besides the point)

u = 10^6; % unit conversion (m to um)
p = 10^12; % unit conversion (s to ps)
% parameters
sigX = 30/u; % radial dimension
sigZ = 1/u; % longitudinal dimension
gm = 100; %energy/gamma
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% array size
zsize = 50;
%parameters
% domain dimensions (um)
zdomain = 5/u;
dz = 2*zdomain/zsize;
z = -zdomain+dz/2:dz:zdomain-dz/2;

% rescaling of position array and parameters
d = gm*z/sigX;
eta = gm*(sigZ/sigX);
% initialisation of array storing integral values
GS = zeros(1,zsize);
% integral kernel
HS = @(x) -1/(4*sqrt(pi)).*(sqrt(pi).*abs(x)+pi*(1-x.^2/2).*exp(x.^2/4).*erfc(abs(x)/2));
% loop evaluates integral for different values of d (s)
for j = 1:1:zsize

    GS(j) = integral(@(x) HS(d(j)-x).*exp(-x.^2/(2*eta^2)),-inf,inf);
        
end

This integrand behaves as follows for d=0, and for the parameter values chosen behaves pretty closely to the kernel, HS.

Running the code, I get NaN values for all d (s), and MATLAB seems to decide very quickly that the integral isn't going to work. I don't understand why but I have the feeling I've missed something obvious. Is there some way in which I've written HS(x) badly?

Thanks!

[Ryan]

Just putting a little comment, sometimes I refer to a variable 's', I mean the array 'z'. Currently the website is not allowing me to edit the post to correct this for whatever reason.

Answer 1

[Paul]

Hi Ryan,

I think that individual terms in the integrand blow up to inf when evaluating in double precision.

u = 10^6; % unit conversion (m to um)
p = 10^12; % unit conversion (s to ps)
% parameters
sigX = 30/u; % radial dimension
sigZ = 1/u; % longitudinal dimension
gm = 100; %energy/gamma
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% array size
zsize = 50;
%parameters
% domain dimensions (um)
zdomain = 5/u;
dz = 2*zdomain/zsize;
z = -zdomain+dz/2:dz:zdomain-dz/2;
% rescaling of position array and parameters
d = gm*z/sigX;
eta = gm*(sigZ/sigX);
% initialisation of array storing integral values
GS = zeros(1,zsize);
% integral kernel
HS = @(x) -1/(4*sqrt(pi)).*(sqrt(pi).*abs(x)+pi*(1-x.^2/2).*exp(x.^2/4).*erfc(abs(x)/2));

Here is the integrand for j = 1

j = 1;
I = @(x) HS(d(j)-x).*exp(-x.^2/(2*eta^2));

Evaluating in double yields a NaN; I think the integrand contains products of very large numbers and very small numbers, and an individual term(s?) blows up before that multiplication when abs(x) is large, which causes a problem when integrating from -inf to inf

I(40)

ans = NaN

Evaluate using vpa arithmetic and we see a reasonable result

double(vpa(I(vpa(40))))

ans = -9.5386e-34

The Symbolic Toolbox can't find a closed form solution (maybe the integrand be simplify-ied first?)

syms x
I = HS(d(j)-x).*exp(-x.^2/(2*eta^2))

I = 

int(I,x,-inf,inf)

ans = 

But does come up with a numerical answer using vpaintegral. I didn't check if that solution is reasonable.

double(vpaintegral(I,x,-inf,inf))

ans = -0.5263

If you want to stay with integral maybe the code can do some up front work to find finite limits of integration where the integrand doesn't become NaN, or maybe a u substitution can be found that keeps the integrand non-NaN and non-inf over the entire real line. Or maybe the integrand can be rewritten/refactored to keep it non-NaN/non-Inf.

[Ryan]

Hi Paul,

Thanks for your comment and suggestions, I've rather awkwardly implimented vpaintegral into my loop and it now yields the result I'm expecting with no NaN values. The code isn't at all pretty, and this can likely be done in a way that's much better but this should be good enough for my purposes as indicated I simply wanted to try and recreate a literature result for comparative purposes and it's now giving me that!

Thank you so very much and Merry Christmas!