numerical problem in integral

Dear all, I have a problem with numerical integral. Basically, the upper bound of this integral should be Inf. However, if I choose the upper bound equaling to Inf, the value of the integral will be NaN +iNaN. Then I change to a small finite number, for an example 3e4, it WORKS.

But if I increase it to 4e4, Matlab returns NaN+iNaN. So I think I have some trouble with the numerical method. In the process of Matlab, one term in the function might be very large and then it causes the trouble. Could anyone help me with this? Thanks

clear all
close all
clc

r0=12.5e-3;            %transmititng coil
r=r0;                   %receiving coil
l=2.5e-3;                 %distance from coil to the conductive medium
c=25e-3;                %thickness of the conductive medium
z=-(l+c);                %position of receiving coil
mu0=4.*pi.*1e-7;
I=1;
sigma4=0;
sigma3_vector=logspace(-2,8,1e2+1);
sigma1=0;
sigma2=0;
f=0.159e9;
omega=2.*pi.*f;

for m=1:1:length(sigma3_vector)
    sigma3=sigma3_vector(m);
fun_A_IV= @(x) besselj(1,x.*r0).*besselj(1,x.*r).*exp(-x.*l).*x...
    .* 2.*sqrt(x.^2+1i.*omega.*mu0.*sigma3).*exp((sqrt(x.^2+1i.*omega.*mu0.*sigma3)+sqrt(x.^2+1i.*omega.*mu0.*sigma4)).*c+sqrt(x.^2+1i.*omega.*mu0.*sigma4).*z) ...
    ./ ((x-sqrt(x.^2+1i.*omega.*mu0.*sigma3)).*(sqrt(x.^2+1i.*omega.*mu0.*sigma3)-sqrt(x.^2+1i.*omega.*mu0.*sigma4))+(x+sqrt(x.^2+1i.*omega.*mu0.*sigma3)).*(sqrt(x.^2+1i.*omega.*mu0.*sigma3)+sqrt(x.^2+1i.*omega.*mu0.*sigma4)).*exp(2*c.*sqrt(x.^2+1i.*omega.*mu0.*sigma3)));

A_IV(m)=integral(fun_A_IV,0,3e4,'RelTol',1e-6,'AbsTol',1e-12,'ArrayValued',true);
end
M=2*pi*r*r0*mu0*A_IV;
figure
semilogx(sigma3_vector,real(M)*1e9,'*')
hold on
semilogx(sigma3_vector,imag(M)*1e9,'o')

[Robert]

Your fun_A_IV returns NaNs for large inputs, especially for larger values of sigma3.

>> fun_A_IV(100,sigma3_vector(70))
ans =
      1.376489474228903e-80 - 6.535524918049053e-80i
>> fun_A_IV(1e6,sigma3_vector(70))
ans =
                NaN +               NaNi

Answer 1

[Walter Roberson]

You are probably getting numeric underflow along with numeric overflow, so you probably have a 0 * inf which would result in NaN.

If you have the symbolic toolbox you might want to switch to int() . That will not necessarily solve the problem, but it would postpone it. (It would not solve it because you are probably going to need to have int() do numeric integration and the symbolic toolbox can underflow too, just not as early.)