I am trying to solve this numerical integration by simpsons method and plot figure. But it returns NAN value.I could not understand why it returns NAN value
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
alpha = 45;
beta = 185;
gamma_e = 116;
G_ei = -4.11;
G_ee = 2.07;
G_sr = -3.30;
G_rs = 0.20;
G_es = 0.77;
G_re = 0.66;
G_se = 7.77;
G_sn = 8.10;
G_esre = G_es*G_sr*G_re;
G_srs = G_sr*G_rs;
G_ese = G_es*G_se;
G_esn = G_es*G_sn;
t_0 = 0.085;
r_e = 0.086;
a = -60;
b = 60;
n = 300;
f = -40:40
w = 2*pi*f; % angular frequency in radian per second
for k = 1:length(w)
L_initial = @(w1) (1-((1i*w1)/alpha))^(-1)* (1-((1i*w1)/beta))^(-1);
L_shift = @(w1) (1-((1i*(w(k)-w1))/alpha))^(-1)* (1-((1i*(w(k)-w1))/beta))^(-1);
L = @(w1) L_initial(w1)*L_shift(w1); low pass filter for w1*(w-w1)
Q_initial = @(w1) (1/(r_e^2))*((1-((1i*w1)/gamma_e))^(2) - (1/(1-G_ei*L_initial(w1)))* (L_initial(w1)*G_ee + (exp(1i*w1*t_0)*(L_initial(w1)^2*G_ese +L_initial(w1).^3*G_esre))/(1-L_initial(w1)^2*G_srs)));
Q_shift = @(w1) (1/(r_e^2))*((1-((1i*(w(k)-w1))/gamma_e))^(2) - (1/(1-G_ei*L_shift(w1)))* (L_shift(w1)*G_ee + (exp(1i*(w(k)-w1)*t_0)*(L_shift(w1)^2*G_ese +L_shift(w1)^3*G_esre))/(1-L_shift(w1)^2*G_srs)));
Q = @(w1) Q_initial(w1)*Q_shift(w1);
p_initial = @(w1) (pi/r_e^4)* (abs((L_initial(w1)^2*G_esn)/((1-L_initial(w1)^2*G_srs)*(1-G_ei*L_initial(w1))))).^2 * abs((atan2((imag(Q_initial(w1))),(real(Q_initial(w1)))))/imag(Q_initial(w1)));
p_shift = @(w1) (pi/r_e^4)* (abs((L_shift(w1)^2*G_esn)/((1-L_shift(w1)^2*G_srs)*(1-G_ei*L_shift(w1))))).^2 * abs((atan2((imag(Q_shift(w1))),(real(Q_shift(w1)))))/imag(Q_shift(w1)));
p = @(w1) p_initial(w1)*p_shift(w1);
I(k) = simprl(p,a,b,n) (simprl is a calling function)
end
plot(f,I)
xlabel('f (frequency in Hz)')
ylabel('I (numerical integral values)')
%%%%%%%%%
function [s] = simprl(p,a,b,n)
h = (b-a)./n;
s1 = 0;
s2=0;
% loop for odd values in the range
for k = 1:n/2;
x = a + h*(2*k-1);
s1 = s1+feval(p,x);
end
% loop for even values in the range
for k = 1:(n/2 - 1);
x = a + h*2*k;
s2 = s2+feval(p,x);
end
s = h*(feval(p,a)+feval(p,b)+4*s1+2*s2)/3;
0 commentaires
Réponse acceptée
Walter Roberson
le 19 Mai 2016
Your code was missing a couple of % and so would not run.
In your line
p_initial = @(w1) (pi/r_e^4)* (abs((L_initial(w1)^2*G_esn)/((1-L_initial(w1)^2*G_srs)*(1-G_ei*L_initial(w1))))).^2 * abs((atan2((imag(Q_initial(w1))),(real(Q_initial(w1)))))/imag(Q_initial(w1)));
when w1 is 0, imag(Q_initial(w1)) is 0 because Q_initial(0) is real, so you have a division by 0.
3 commentaires
Walter Roberson
le 20 Mai 2016
You have
f = -40:40
w = 2*pi*f;
so f passes exactly through 0, and at the point it is exactly 0, w will be exactly 0, which leads to that problem about division by 0.
So if you add a tiny amount to f so that it is never exactly 0, then I think it might do the trick.
Example:
f = (-40:40) + 1/10000;
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Food Sciences dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!