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?
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
Namira
el 28 de Sept. de 2016
Comentada: Namira
el 28 de Sept. de 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;
Respuesta aceptada
Walter Roberson
el 28 de Sept. de 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.
3 comentarios
Walter Roberson
el 28 de Sept. de 2016
P0vals = [6 -6];
for P0idx = 1 : length(P0vals)
P0 = P0vals(P0idx);
...
I(k, P0idx) = ...
end
Más respuestas (0)
Ver también
Categorías
Más información sobre Frequency Transformations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!