MATLAB: Infinite or Not-a-Number function value encountered.

Question

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

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.