Infinite or Not-a-Number function value encountered.

Good evening,

I would like to plot (t,f(t)), with

 f(t)=integral_0^ (2*pi*1.0e+12)*75.2933 y(t,w) dw

For that, i used 'quad' to compute my integral, but there is a problem when matlab run: Warning: Infinite or Not-a-Number function value encountered.

Thanks i advantage for you help.

function y=toltallA(w,T)
%%%Parameter
%%Begin parameter
L=2.9e-6;
v= 21.6e+3;
a=0.142e-9;
N_0=6.022e+23;
A= 3*sqrt(3)*a^2*N_0/4;%
w_max=(2*pi*1.0e+12)*75.2933;
delta_M=3;
B_N=3.85e-25; %%
B_U=7.7e-25;%%
alpha=3; %%
theta=1000;
c_d=1.0e-6;
h=6.626e-34;
hbar= h/(2*pi);
kB=1.38*1.0e-23;
coef=hbar^2/(4*pi*kB*v^2);
%end of my parameter

%%%L'unique fonction qui depend de ''w''. J'integre par raport a w et T est un parametre

tau=1./(v/L+ (A*delta_M^2*2*pi*c_d)/(2*pi*v^2*w_max^2)*w.^4+(B_N+B _U*exp(-theta/(alpha*T)))*T^3*w.^2);%%
y= coef*v^2/T^2*tau.*w.^3.*exp(hbar*w./(kB*T))./(exp(hbar*w./(kB*T))-1).^2;

%%%MAin Program j'utilise quad pour l'aproximation de l'integral et loop en T.
T=1:10:1000;
QlA=zeros(size(T));
for k=1:length(T)
QlA(k)= quad(@(Z)toltallA(Z,T(k)),1.e-20,(2*pi*1.0e+12)*75);
end

Answer 1

[Roger Stafford]

As I stated in your similar query in the matlab newsgroup, your trouble occurs because the exponent, hbar*w./(kB*T), becomes excessively large causing matlab's 'exp' function to overflow to infinity and then the ratio in y to produce a NaN. Besides the remedy I stated there you can also use

 1/4*csch(hbar*w./(2*kB*T)).^2

instead of

 exp(hbar*w./(kB*T))./(exp(hbar*w./(kB*T))-1).^2

since these are also equivalent.

You may also run into problems at the lower end of w values. When w is very small the above quantity approaches infinity, though it is multiplied in the y expression by w^3 which more than compensates. However, to maintain accuracy and/or avoid another NaN you may have to approximate the expression for y using a simple Taylor approximation of csch(x)^2 with w in the neighborhood of zero.

Answer 2

[julien]

Thank you very much for the idea. I tried to remplace f=w.^2.*exp(hbar*w./(kB*T))./(exp(hbar*w./(kB*T))-1).^2 by a funtion like: function f=fcond(w,hbar,kB,T) f(w<1.e-2*kB*T/hbar)=1; %%%one is the approximation of xx on zero.

f(1.e-2*kB*T/hbar<=w)=1/4*w(1.e-2*kB*T/hbar<=w).^2.*csch(hbar*w(1.e-2*kB*T/hbar<=w)/(2*kB*T)).^2;

My idea is : if hbar*w./(kB*T) <=1.e-2 we take f=1; otherwise your approximation for hbar*w./(kB*T) =>1.e-2

What do you think about my function; and how to inert this function in the main program. Many thanks.

[julien]

Thank you very much for the idea. I tried to remplace f=w.^2.*exp(hbar*w./(kB*T))./(exp(hbar*w./(kB*T))-1).^2 by a funtion like: function f=fcond(w,hbar,kB,T) f(w<1.e-2*kB*T/hbar)=1; %%%one is the approximation of xx on zero.

f(1.e-2*kB*T/hbar<=w)=1/4*w(1.e-2*kB*T/hbar<=w).^2.*csch(hbar*w(1.e-2*kB*T/hbar<=w)/(2*kB*T)).^2;

My idea is : if hbar*w./(kB*T) <=1.e-2 we take f=1; otherwise your approximation for hbar*w./(kB*T) =>1.e-2

What do you think about my function; and how to inert this function in the main program. Many thanks.

[Roger Stafford]

It might be simpler to just integrate over two separate w-intervals with one very near zero and the other all the remaining w-values, instead of writing a conditional function. Either way ought to work, however.

The function you are integrating is of the form

 y = k1./(k2+k3*w.^2+k4*w.^4).*w.^3.*csch(k5*w).^2 =
 k1./(k2+k3*w.^2+k4*w.^4).*w.^3./(k5*w+(k5*w).^3/6+(k5*w).^5/120+...).^2 =
 k1./(k2+k3*w.^2+k4*w.^4).*w./(k5+k5^3*w.^2/6+k5^5*w.^4/120+...).^2

where the second denominator is the Taylor series for sinh(k5*w). For the first w-interval you can approximate this by having only a few of these Taylor series terms. The more terms you include, the larger you can make this interval and therefore the less error you have in the y computation for the second interval where you have to use csch(k5*w) and are threatened by a zero divided by zero situation when w is near zero. You should choose some suitable compromise.

For numerical accuracy you should probably not be using the 'quad' function to do the integration, but rather 'quadgk', 'quadl', or better still 'integral' since these can better handle awkward integrands such as you have here.