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

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;

Answer 1

[Walter Roberson]

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.

[Namira]

Thanks. How can I solve it?

[Walter Roberson]

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;

[Namira]

Thanks. @ Walter Roberson