# I am trying to solve this integration by simpsons method and plot figure. It's give figure taking only the value of P0= -6 from for loop. How can I solve it?

3 views (last 30 days)
Namira on 28 Sep 2016
Commented: Namira on 28 Sep 2016
alpha = 45;
beta = 185;
gamma_e = 116;
% Gain values
G_ei = -18.96;
G_ee = 18.52;
G_sr = -0.26;
G_rs = 16.92;
G_es = 2.55;
G_re = 4.67;
G_se = 0.73;
G_sn = 2.78;
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; % corticothalamic loop delay in second
r_e = 0.086; % Excitatory axon range in metre
f = linspace(-40,40,500); % f = frequency in Hz
w = 2*pi*f; % angular frequency in radian per second
delt_P = 0.5;
L=zeros(1,500);
Q=repmat(L,1);
P=repmat(L,1);
%%%%%%%%%%%%%integration %%%%%%%%%%%%
a = -80*pi;
b = 80*pi;
n=500;
I=repmat(L,1);
P_initial = repmat(L,1);
P_shift = repmat(L,1);
p = repmat(L,1);
for k = 1:length(w)
for P0 = [6 -6]
L_initial = @(w1) (1-((1i*w1)/alpha))^(-1)*(1-((1i*w1)/beta))^(-1);
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)));
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)));
G = 150*exp(- (f - P0).^2./(2*(delt_P).^2));
P2 = @(w1) G(k) + P_initial(w1);
L_shift = @(w1) (1-((1i*(w(k)-w1))/alpha))^(-1)* (1-((1i*(w(k)-w1))/beta))^(-1);
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)));
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) P2(w1)*P_shift(w1); % Power spectrum formula for P(w1)*p(w-w1)
I(k) = simprl(p,a,b,n);
end
end
figure(1)
plot(f,I,'r--')
figure(2)
plot(f,G,'k')
%%%%% Simprl %%%
function [s] = simprl(p,a,b,n) % f is the function to be integrated; a = initial value of the interval; b = final value of the interval; n = No. of subintervals
% The function implements the Simpson's Rule
h = (b-a)./n;
s1 = 0; % The variable s1 is initialised to 0.
s2=0; % The variable s2 is initialised to 0.
% loop for odd values in the range
for k = 1:n/2; % The index variable k starts at 1, then increases in steps of 1 until it reaches n/2.
x = a + h*(2*k-1);
s1 = s1+feval(p,x); % Each time through the loop the value of feval(f,x) is added to s1.
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
% Final result of integration where odd values are multiplied by 4 and even values are multiplied by 2
s = h*(feval(p,a)+feval(p,b)+4*s1+2*s2)/3;
Walter Roberson on 28 Sep 2016

Walter Roberson on 28 Sep 2016
You have the structure
for k = something
for P0 = something
...
I(k) = simprl(p,a,b,n);
end
end
Each time you iterate through the different P0's, you store only one value, the value for the active k. The effect is as if you had only done the calculation for the last P0.
You should be indexing your I according to k and also according to how many P0 you have done.
I(k, P0idx) = ...
where P0idx would be something you had set at 1 for the first P0 index, at 2 for the next P0 index, and so on.
Namira on 28 Sep 2016
Thanks. It helps me a lot.