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

2 views (last 30 days)
Namira on 19 May 2016
Edited: Namira on 22 May 2016
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;

Walter Roberson on 19 May 2016
Your code was missing a couple of % and so would not run.
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.
Namira on 22 May 2016
Thanks. @ Walter Roberson