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 次查看(过去 30 天)
显示 更早的评论
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
2016-9-28
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 个评论
Walter Roberson
2016-9-28
P0vals = [6 -6];
for P0idx = 1 : length(P0vals)
P0 = P0vals(P0idx);
...
I(k, P0idx) = ...
end
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Frequency Transformations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!