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
采纳的回答
Roger Stafford
2013-4-3
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.
更多回答（1 个）
julien
2013-4-7
0 个投票
2 个评论
julien
2013-4-7
Roger Stafford
2013-4-7
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.
类别
在 帮助中心 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
