通过辛普森方法数值积分

Question

我想通过辛普森方法来解决这种积分，并将结果绘制在一个图中。该图仅取自for循环的P0 = -6的值。当我设置我（K，P）时会给出错误：

Attempted to access P0(0); index must be a positive integer or logical

我该如何解决？

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') 

Answer 1

目前只为你将它们保存在I(k)使用P0 = -6结果。首先你将结果保存为P0 = 6，稍后覆盖它并保存另一个。 P0 = 6的结果既不使用也不保存。如果您按照以下方式编写代码，则需要澄清。

for k = 1:length(w) 

    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)));   


    for P0 = [6 -6] 
     G = 150*exp(- (f - P0).^2./(2*(delt_P).^2)); 
     P2 = @(w1) G(k) + P_initial(w1); 
     p = @(w1) P2(w1)*P_shift(w1);   
     I(k) = simprl(p,a,b,n); 
    end 
end 

您不能访问I(k,P)为I是一个向量不是一个矩阵。但是这会给你索引超出矩阵的尺寸。您可以将P0 = -6的结果保存在一个变量中，P0 = 6保存在另一个变量中，因为您的代码中的结果不会相互依赖。